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Abstract 


This note is part of a continuing study of future problem, areas in 
structural dynamics of space vehicles, conducted by the author for the 
Jet Propulsion Laboratory. 

The motivation for this particular piece of work is the conviction 
that future space vehicles will be relatively large and flexible, and that 
active control will be necessary to maintain geometrical configuration. 
While the stresses and strains in these new space vehicles are not expected 
to be excessively large, their cumulative effects will cause significant 
geometrical nonlinearities to appear in the equations of motion, in addition 
to the nonlinearities caused by material properties. Since the only ef- 
fective tool for the analysis of such large complex structures is the digital 
computer, it will be necessary to gain a better understanding of the non- 
linear ordinary difference equations which result from the time 
discretization of the semi -discrete equations of motion for such structures. 


1 , Intr odu cti on 


Equations of the type; 


X , = f (x , n) 
— n-fl — n ’ 


X = c 
— o — 




X = c 
— o — 


are Imown as nonlinear ordinary difference equations or point mappings. 
Equation (l, l) is known as an explicit nonlinear difference equation, while 
Eq. (1.2) is known as an implicit nonlinear difference equation. 

If in Eq. (1 . l) 

0 

f a) = A(n)x^+ g(n) (1.3) 

then (1.1) becomes 


^n+1^ if"') 


X = c 
o — 


Similarly, if in (1.2) 


l(^En'-n+l’“^ ^ A{n)x^+ B(n)x^_^^ + g(n) 


Then (1.2) becomes 


in+r 


X = c 


{ 1 - 6 ) 


Equation (1.4) is Icnown as a linear explicit difference equation, while (1.6) 


" 3 - 


is known as an implicit linear difference equation. 

Since Eq. (1.6) can be rewritten as: 

5 „+r C(n)x^+ h(n) 

X = c 

^ - . (1.7) 

C(n) = [I-B(n)]"^A(n) 

b(n) = [I-B(n)]"^g{n) 

Thus, there is no difference, in theory, between ejqplicit linear difference 
equations and implicit linear difference equations. Unfortunately the same 
is not true, in general, for nonlinear difference equations. 

Difference equations arise in a variety of scientific and engineering 
disciplines, for example; 

(a) In biology; population genetics and dynamics are described by 
nonlinear difference equations, 

(b) In control theory; sampled data control system are described 
by either linear or nonlinear difference equations. 

tl 

(c) In numerical analysis; in order to solve a differential equation 
on a digital computer, the independent variable must be discretized 
and the differential equation becomes a difference equation. In 
particular, nonlinear differential equations become nonlinear dif- 
ference equations. 

It is to this last class of problem that this note is addressed. 

2, Existence and Uniqueness of a Solution of the Initial Value Problem 
a) Explicit Nonlinear Difference Equations 

Theorem 1 Given the explicit nonlinear difference equation 


-4- 



X = c 
— o — 




( 2 . 1 ) 


If (i) Vx, £(x) is continuous in x, therefore l]f(x)|l <co, Vj^H^oo 

(ii) ||cll<oo. 

Then there exists a uriique solution of the initial value problem {2, l) 

Proof 

Since lyi <oo 


llxj 11^ 


Therefore there exists a solution to Eq, (2, l), satisfying the initial data. 

Since the process of generating the solution is explicit, then there exists 
one and only one solution of {2. 1) satisfying the initial data, therefore the 
solution of (2. 1) is unique. It will be noted that for explicit nonlinear differ- 
ence equations, the question of existence and uniqueness of a solution is S 

trivially answered in comparison with the same question for nonlinear dif- ^ 

ferential equations. 

b) Implicit Nonlinear Difference Equations 

Let us consider nov’ the implicit nonlinear difference equation 

^ . 


X = c 

— o — 


( 2 . 2 ) 
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In the general case, we can say relatively little about the existence o£ a 
solution to Eq, (2.2). The implicit function theorem guarantees, under weak 
restrictions on f (x ,x ,), that there exists a unique local solution of 
(2.2) provided 1 jclj is sufficiently small. In some special cases Eq. (2.2) 
may be inverted so that it is described by an explicit equation. 




= c 



(2.3) 


In the case of most practical importance, Eq. (2.2) has the 
structure 


X X + ef ,(x ,x ,t) 
— n+1 — n — 1'— n— n+l'' 


X = c 


(2.4) 


where |e 1 is frequently a small quantity. 

Before proving the existence of a unique solution of Eq. (2.4) we 
will establish the following theorem. 

Theorem 2 Given the implicit equation 


x = g(x) (2.5) 

and the iterative procedure 

= g{^) n = 0,1,2 .. . . (2.6) 

Then if ^(x) satisfies the following conditions 

(i) llg(3£) - g(y)!l llx-vj| £orVx,yeS S: lj 2 i-x° || £ p 
, with 0 ^ \ ^ 1 

(ii) There exists an 5 ijg(x°)-x°||£ (l-\)p 
then V iterates x^ satisfy the following conditions 


(2.7) 

( 2 . 8 ) 
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(i) lkn'25oll 


(ii) lim llx 11 = Of where a - g(oi) 



(2.9) 

(iii) Of is the only root of Eq. (2.5) in j 

J 



Proof 




Since x - g(x ) 

— n+1 ~'_n' 



(2.10) 




(2,11) 

•■• lhn+rSil>= 'fefen)-sK-i)il 



(2. 12) 

£ X 1 lx - X 11 
' t-n — n-i “ 



(2.13) 


Now llxj-x^ll s(l-\)p <p (2.14) 

X ^ € S 

•'• !fe“-iil + n^Er^ol 1 

^[X(l-X) +(1-X)]p =(l-x^)p<p (E.16) 

Suppose that x , x, ... x £ S 
— o —1 “n 



(2.17) 


(2. 18) 

fc„+rioll^fen+r5n"+fen-in-lll+ ' • • 

(2.19) 

+ . . . X+l)(l-X)p 

(2.20) 

s(l-x“'^^)p<p 

(2.21) 

5n+l « S 



I 4 
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fen+k-2j^fen+h"^n+h-lll+ • ‘ • • + fen+l'^J 

^/,n+h-l .n+h-2 , .iiw^ .. 

^(A. +X + ... X )(1-X)p 

^ A P 


■\^S- '^-n+k-Si"=0 

.*. the sequence (x^l is a Cauchy sequence and converges uniforinly, 

L,im 3£ , = O’ = L/im g(x ) = g (Litn x ) = gfo) 

n“*co~n+l — n“+oo — n' ^'‘n*»oo n'' 


( 2 . 22 ) 

(2.23) 

(2.24) 

(2.25) 

( 2 . 26 ) 
(2.27) 


the sequence (x^l converges uniformly to a limit, a 6S, which is a 
solution of Eq. (2.5). 


Uni city 


If Of and p are solutions of Eq, 
Then 

(2.5) which both belong to the 

1« 

II 

o es 

(2.28) 

fi = g(£) 

P €S 

(2.29) 

11“ “ fill 

1 = y^) - g(p)ii^xiia - pii 

(2.30) 

•: ilH - £ll 

(l-X)^O £^P . 

(2.31) 


Thus, under the hypothesis of Theorem 2 there exists a unique solution of 
Eq. (2.5). 

■ Returning to the question of the existence and uniqueness of a 
solution of the initial value problem for an implicit nonlinear difference 
equation, we have Theorem 3. 



- 8 - 


Theorem 3 Given the implicit nonlinear difference equation 


— n+1 — n — '-^n'— n+i' 


X 

— o 


= c 


(2.3Z) 


If 

(i) l(HEn'-n+l^ is continuous in ^ and 

(ii) continous first partial derivatives with respect 

to X 
— n+x 

(iii) |e| is sufficiently small. 

Then there exists a unique solution to Eq, (2.32) on some finite interval 
0 £ n < N. 

Proof Let x = = Sn ^ 5n+l^ "" 

With hypothesis (i) and (ii) g(x) satisfies 


Ikfe) ■ g(z)ll= k I liltin’ S) ‘Ifen’^MI 

II J( 5)11- lli-y|l 

£X lix-jrll 


where 

J(S) =1- Sitin’ S) lx =? 

^ = ofx + ( 1 -or)y 0 < or< 1 


(2.34) 


r 

(2.35) 


I 


For 

^ te " 2 eJI ^ llx"Enll^p 


(2.36) 


We can always choose &\ sufficiently small so that 


-9- 


'■= unu(i)ii<i 

If we choose = x as oui- initial iterate, then 

= II ®i(^n^^n^ll 


(E.37) 


(2.38) 


Since f(x,y) is continuous in x and y we can always choose je[ suf- 
ficiently small so that 

"^°II^ (1“X)p (2.39) 


Thus, given the hypothesis (i) and (ii) we can always choose | ej sufficiently 
small so that the conditions of Theorem 2 are satisfied. Thus given an x^, 
there exists a unique solution satisfying 


^+1 + ( 2 - 40 ) 

Thus, sh.rting with x = c and lej fixed and sufficiently small, 
there exists a unique x , satisfying Eq, {2.32). If, using x and the 

■i 

same value of e, conditions (2.37) and (2.39) are satisfied, then there 
exists a unique x^ satisfying Eq. (2.32). Proceeding in this way, we 
check at each step to see if conditions (2.37) and (2.39) are satisfied. If 
they are satisfied at each step, the solution can be continued indefinitely 
into the future. If they are not satisfied after a finite number of steps 
the solution may cease to exist or go to infinity. Thus, given condition 
(i), (ii) and (iii) there exists a unique solution to Eq. (2.32), at least 
on some finite interval 0 < n <N. 

Theorems 1 and 3 deal with autonomous equations that is, equa- 
tions which do not contain n explicitly. The hypothesis of Theorems 1 
and 2 can be relaxed to include explicit dependence on n, in addition to 
domain dependent continuity properties. 


- 10 - 


3. Properties of Linear Diffei’ence Equations 


(a) Diflereaace Equations with Constant Coefficients 
Consider the linear difference equations 

(3,1) 

Ho “ £ l-A-1 oj 

where A is a constant matrix witla |1 a|| = a “^co, Hf(n)H'^oo, The solution 
of Eq, (3. l) is easily formed hy elementai*y metliods 


Hn+l" ^n 


» « 


X. Ac t f ^ 

— 1 — — o 

2 

X, = Ax +f, = Ac+Af +f^ 

— i —1 —1 — — o —1 

2 

x^ = + f_2 = A^ c + ^ A^ 

i=o ^ 

n 


X 


■n+1 


= 


y 

i=o 




A^"^f . 
—1 


(3.2) 

(3.3) 

(3.4) 


(3.5) 


Alternatively we can write tliis solution in terms of tlie principal matrix 

solution X , where 
n 


thus 


X ^ = AX ; X = I 
n+1 n o 


= A 


X^ =A‘ 


(3.6) 

(3.7) 

(3.8) 


X^ ^ A 

II 


n 


(3.9) 


-u- 


Thiis, the solution to Eq. (3.1) can be written 

X ,.“X ,,c+) X .f. 

— n+1 n+1— .Ll n-i — i 

1=0 

We note that X . = 
n-i 


n 


-n+1 


= X 


n+l-+Z 


X f. 

n+1 x+1 — 1 


1=0 


(3.10) 

(3.11) 
(3.: 2) 


("b) D^f^<srence Equations with Variable Coefficients 
Consider the linear difference equation 

X . - = A(n)x + f 
—n+1 ^ n — n 

> 

x^ = c [a| 

where A(n) is a step dependent matrix, with |jA(n)|j< oo, Vn and 
j|£(n)ll <oo, Vn. 


(3.13) 


The solution of Eq. (3. 13) is also easily formed by elementary 
methods. 

x^=A(0)c+f^ (3.14) 

£2 = A(l)x^ +f_^ = A(l)A(0)c + A(l)f^+£^ (3.15) 

x^ = A(2 )x2 +I 2 = A(2)A(l)A(0)c + A(2)A(l)f^+A(2)f ^ +f 2 (3.16) 

r n i n 


-n+1 ^ A(n)A(n-l) 


A(0) 


£+1 If A(h)"^f. 

i=o n=o 


(3.17) 


Alternatively we can write tliis solution in terms of the principal matrix 

solution X , where 
n 

X^ = A(n)X , X =I 
n+1 ' ^ n ' o 


(3,18) 


~ 12 - 


a 


Thus 

“ A(0) (3*19) 

X2=A(1)A(0) (3.20) 


• n-1 

X = A(n-l)A{n-2) . . .A(0) A(i) (3.21) 

i=o 


Thus, the solution to Eq. (3.13) can be written 

n 


-n+1 “ ^n+1' 


^ -"I 


X , , X."}. f . 
n+1 1+1 “1 


1=0 


(3.22) 


In general we can say very little about the structure of the solution in the 
case of variable coefficients. There is, however, one special case, the 
case of a difference equation with periodic coefficients. 


(c) Difference Equation with Periodic Coefficients 
Theorem 4 Consider the homogeneous difference equation 


X , = A(n) X 

n+i '■ ‘ n 


X = I 
o 




where 

A(n+N) = A{n) |A(n) j 0 
||A(n)lj<o' Vn 





The principal matrix solution X^ has the form 

Xjj =Q(n)c“ 


(3.23) 


(3.24) 


(3.25) 


-13- 


wiiere Q(n+ISf) ?=Q(n) is a periodic matrix and C is a non- singular con- 
stant matrix. 

Proof Prom (3.23) 


X = I . , 
o 

(3.26) 

Xj^ = A(0) 

(3.27) 

X^ = A{1)A(0) 

(3.28) 


N-1 

^n=]l 

1=0 

(3.29) 

We note that Xj,^ = Jj A(i) i s non-singular since A(i) 

i=o 

Vi. 

is non- singular 

^+1 =n 

i~o 1=0 

(3.30) 

N+1 1. 

^+2=Jl 

1=0 1=0 

(3.31) 

Ntlc-1 krl 

^N+k”fl ~n 

i=o i=o 

(3.32) 

Similarly 



(3.33) 


Since Xj^ 



is non- singular, we may write 
K 

= G , C - a constant matrix 


(3.34) 


Consider the matrix 


-14- 


Q(^) = 


-n 


Since 


= I , Q(0) = I 


Thus 


Q(N+n) 


-(N+n) 

•^N+n^ 


= X X-^C'^C~^ 
n N 


But 


-N 

i 


.*.Q(N+n) = =Q(n) 


Q (n) is a periodic mairix with, period N. 


(3.35) 


(3.36) 


(3.37) 

(3.38) 


(3.39) 


Hence 

X^=Q(n)C^ (3.40) 

4 

Thus the complete structure of X^ is known for Vn if Xj^ is known for 
0 ^ h £ N. 

We note in passing that if A is a constant matrix, then A is a 

i 

periodic matrix of period N = 1, hence, difference equations with constant 
coefficients are a special case of difference equations with periodic co- 
efficients and in this case the matrix C = A and Eq. (3.40) becomes 

X^ = a” (3.41) 

4. Stability of Difference Equations 
Definition Liapunov Stability (E.S.) 


Given the difference equation 


-15- 

5a+l = £<2n’ (^-1) 

where 

l(£^£)^i (4-2) 

The ecjihlibrimn solution ^ = 0 is said to be Itiapunov stable if 
given any 6>0, there esdsts an e>0, such that if jjx ||<e, then 
Ijx^ll < 6 for all n > 0 . 

Liapunov Asymptotic Stability (ij.A.S. ) 

The equilibrium solution x = 0 is said to be Liapunov asymptotically 
stable if (a) it is Liapunov stable and (b) 11 ^^^11"* 0 as n -• oo. 

(a) Stability of Linear Difference Equations 

(i) Linear Difference Equations with Constant Coefficients 

Theorem 5 Given the difference equation 

X - Ax (4, 3) 

-n+l — n 

A “ a constant matrix 

(i) If A is non-defective (i.e. has a full complement of ordinary 
eigenvector) necessary and sufficient conditions for Liapunov 
stability are that the eigenvalues of A should be less than or 
equal to imity in mod'olus. 

(ii) If A is defective (i. e. does not have a full complement of 
ordinary eigenvectors) necessary and sufficient conditions for 
Liapunov stability are that the eigenvalues of A should be less 
than unity in modulus. 

Proof 

(i) If A is simple, ie. non-defective there exists a similarity 



-16- 


matrix 

T 3 T"^AT = A 

where A is a diagonal matrisc. 

A has the representation 

A = TAt"^ 

As previously shown, the principal matrix solution of (4.3) isi 

= A^ = (TAT”^)^ = TA^ = TA^T'^ (4.4) 

Sufficien cy 

If 1X.(A)Hi, then \f"(A) remains bounded as n co 

X is bounded and remains bounded as n "*oo 
n 

A l|xJ^M<oo , Vn (4.5) 


From (4. 3) 



X = X X 
— n n— o 

(4. 6) 

if 

li 2£oil ^ ^ 

(4.7) 

» 

II 2„ll ^ llxj e 

(4. 8) 


^ Me 

(4.9) 


e ^ 5/M 



11x^11 ^ 6 Vn 

(4. 10) 

(4.3) 

is Liiapunov stable at x = 0. 


Necessity 



If 1X^(A)|>1 for some i, then X^(A) cannot remain bounded as 


-17- 


iX'^CD, Hence X . cannot remain bounded ag n"*oo, 
n 

(ii) if A is defective, it cannot be reduced to diagonal form, 
however, there is a similarity matrix T 3 


T"-^AT = 


a. 


J . 


OJi 


= J 


(4. 11) 


where the are Jordan blocks associated with the eigenvalues X.{A) 


ie{l,h) A has the representation 
A - TJT"^ 


(4.12) 


As previously shown, the principal matrix solution of (4.3) is: 


= A^ = (TJT~^‘)”' =: TJ^T~^ 


where 





“'h 


(4.13) 


(4. 14) 


and 


-18- 


and 


J . 

oa 





n(n-l) 

i 



nX^"^ 

1 

1 


1 


( 4 . 15 ) 


L . 

Sufficiency If 1 Xj(A) | < 1, then J^, remains bounded and tends to 

^ °^i 

zero as n -^oo. X is bounded and tends to zero as n*^oo. 

n 


^ oo Vn 
and Lim||xJ|-»0 


(4.16) 


n-^oo 


From which we immediately deduce that if |X^{A) j < 1 Vi, the system (4. 3) 
is not only Liapunov stable, but is asymptotically stable. 

Necessity If j\.(A)j^ 1 for some i, then cannot remain bo^mded 

“i 

as n*^co, hence X cannot remain bounded as n-^oo. 

n 

Alternatively use can be made of Liapunov's Theorem. 

Theorem 6 (Liapunov) Given the difference equation 


X , = Ax 
— n+1 — n 


A - a constant matrix 

Then (4. 17) is Liapunov asymptotically stable at x = 0 iff there exists 
a symmetric positive definite matrix P such that 


(4.17) 


A-^PA - P = - Q 


(4.18) 


•19“ 


Proof Stifficfency 

Suppose that there exists a matrix P satisfy (4. 18) let 


P X 

n —XI — n 


(4.19) 


Since P is symmetric and positive definite is positive definite 


T „ 

v , T = X P X 
Vn+1 -n+i -n 


(4.20) 


Using (4. 17) 




n*fl 


—n' ■ — n' 

T T 

= x„A PAx„ 

— n — n 


= ’^n+r'^n = (A^PA-P)x^ 


(4.21) 

(4.22) 

(4.23) 


Using (4. 13) 




(4.24) 


Thus 


V <V < V <V , . . . < V 
n+1 n n-1 n-2 o 


(4.25) 


Since vanishes only at the origin 


V 0 as n CO 
n 


(4.26) 


If we define 


II X 11 =Jy 

— n"n n 


(4.27) 


we see that 


X < jx 

r“\ 


■— n 


° P 


(4.28) 
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If 

Ilx^|ls£sx6 (4.29) 

then 

!lxjlp<6 vn (4.30) 

and 

as n-*oo (4.3l) 

Eq. {4. 17) is Liapxmov asymptotically stable at the origin. 

Necessity 

Let A be a stability matrix e. |\^{A) I < 1 Vi. Let P satisfy 
(4. 18) i. e. 

A*^PA - P = -Q (4.32 ) 


We wish to show that P is syinmetric and positive definite. 

T 

If we premultiply (4.20) by A and post multiply by A, then 
2T 2 

by A and A , etc., we obtain 

"\ 

A*^PA - P = - Q 

2T 2 T T 

A PA -A PA ^ - A QA 

. 2 T .^.2 . . 2 T ^,2 

(4.33) 



A^ PA^-A^“^ 


PA 


n-1 


= - A 


n-i 


QA 


n-1 




Adding, we obtain 
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n-1 

(A^)'^PA^-P = "5^ (A^)^QA^ (4.34) 

i=o 

Since A is a stability matrix A^-*[0] as u oo 
oo 

P=:^ (A^)'^QA^ (4.35) 

i=o 


We note that; 


i that; T 

/oo \ ^ CD 

(1) p'^iy (a^)'^qam =y (a^)'^q'^a^ 
\uo I i=o 


But Q-" = Q 

= P 


03 


(2) x*^I^=y (A^x)'^Q(A^x) 


1=0 


But Q is positive definite 


(A^x)'^Q(A^x) > 0 


(4.36) 

(4.37) 

(4.38) 

(4.39) 


provided 


A^x ^ 0, 


© If 

1a| ^ 0 Ia^It^O x?tO (4.40) 

x*^P^>0 xj^O (4.41) 

Thus, if A is a stability matrix there exists a P, symmetric 
and positive definite such that Eq. (4.32) is satisfied. 

@ Note If jA| = 0, it appears that (4.41) is not satisfied. How- 
ever, if 1 a| = 0, then A has one or more zero eigenvalues, and the dis- 
placements in these modes vanish after one step, thus the problem is 
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really oiie in (N-h.) dimensions, where h is the multiplicity of the zero 
eigenvalue. Thus, if in (4.29) x 6 IL, the range space of A, P is 
positive definite. 

(ii) Linear Difference Equations with Periodic Coefficients 
Theorem 7 Given the difference equation 


2n+l " 


— n 


A(n+N) = A(n) 


(4.42) 


(i) If the principal matrix solution is simple, necessary and sufficient 
conditions for Liapunov stability are that the eigenvalues of should be 
less than or equal to unity in modulus, 

(ii) If the principal matrix solution is defective , necessary and 
sufficient conditions for stability are that the eigenvalues of X^^ should be 
less than unity in modulus. 

Proof As previously shown, the solution of (4.42) with initial data 
X = c is given by 


X = X c 

— n n— 


(4.43 ) 


where 

= Q(n)C^ (4.44) 

Q(n+N) = 'Q(n) is a periodic matrix and G - X^^^ is a -rionstant matrix. 

If Xp^ is simple, there exists a similarity matrix T 3 

T ^Xp^T = A a diagonal matrix (4.45) 

X^ therefore has the representation 
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= tat“^ 

(4.46 ) 

n 

r. tA^t"^ 

(4.47) 


If clearly and hence renaains hounded as n-+oo, 

therefore (4.42 ) is Liapunov stable at x = 0. 

The remainder of the proof closely follow that of Theorem 5 and will 
not he repeated here. 

Theorem 8 Theorem 6 can be generalized to the case of linear differ- 
ence equations with periodic coefficients. 

Given the difference equation 


X , = A(n)x 

— n-fl ' n 


A(n-fN) = A(n) , jA(n) [ ^ 0 |jA(n)H<oo Vn 


(4.48) 


Then (4.46) is Liapunov asymptotically stable at x = 0 iff there 
exists a symmetric positive definite periodic matrix P(k) such that 


i) P(k+N) = P(k) = p"^(k) positive definite 

ii) A'^(k)P(k+l)A(k) - P(k) = -Q(k) 

T 

iii) Q (k) = Q(k) = Q(k+N) positive definite 


(4.49) 


Proof Sufficiency 

Suppose that there exists a matrix P(k) satisfying (4.49). Let 

(4. SO) 


Since P(n) is symmetric and positive definite for all n, is positive 
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defisixte . 

^n+l " 

Using (4.48 ) 

^n+l “ 

AV^ = ■^n+1 -■^n = 


Using (4 .490 ^i) 


AV^ = -xjQ(n)x^ <0 
n — n ' n. 


Since vanishes only at the origin 


V “*0 as n "♦ CO 
n 


.*. Equation (4,48) is Liapunov asynaptotically stable at the origin. 


(4.51) 


(4.52) 


(4.53 


(4.54) 


(4.55 ) 


Necessity 


Let A he a stability matrix so that V soluLions of (4.48 ) tend to 
zero as t "* co. 

Let P{k) satisfy (4.49) e, 

A'^(k)P(k+l)A(k) - P(k) = -Q(k) (4.56) 

similarly 

A'^(k+l)P(kt2)A(k+i) - P(k+1) “ -Q(k+1) (4.57) 

T 

If (4.57) is premultiplied by A (k) and post multiplied by A(k) we obtain 

(A{k+i)A(k))^p(k+2)(A(k+l)A(k)) - A(k)'^P(k+l)A(k) = - A(k)^Q(k+l)A(k) 

(4.58) 
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Similarly 

A'^(k+2)P(k+3)A(k+2) - P(k+2) - - Q(k+2) (4.59) 

T 

If ( 4.59 ) is premnltiplied by (A(k-H)A(k)) and postmultiplied by 
A(k+i)A(k), we obtain 


( A(k+2)A(k-f-l)A(k)? P{k+3)(A{k+2)A{k+i)A(k) 

- (A{k+l)A(k)'^P{k+2)(A(k+l)A(k)) 

- - (A(k+l}A{k)^Q(k+2)A(k+l)A(k)) 

(4,60) 


Repeating the procedure n times gives 



n^E 

P(k4n-l)^j A{k+i)| 
i=o 



(4.61) 

If these n equations are added, we find that just as in (4.33 ) we obtain can- 
cellation in pairs and finally we have: 


A{k+i)| P(k+n) j jf A(k+i) 


1=0 


/ 


\i=o 

/j.l 


- P(k) 


= -1; 

j=o 


A(k+i) rQ(lrtj)r|f A(k+i)| (4.62) 


1=0 


1=0 


I 


Now 


f 


1=0 


A{K+i)= ^(k+n,k) 


(4.63) 
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where 

$(m,k) = 

satisfies the equation 

^(n+l,k} = A(n)'^(n,k) 

#(k,k)=I 

Since V solutions of (4.48) tend to sero as n^oo 
^{m, k) “*0 as m do 
Thus as n-»oo Eq. (4.62) becomes 
oo 

P(k) «(j,k)'^Q(3)9{j.i^) 

j=k 

OO 

i) P(k)'^=y^ $(j,k)'^Q(j)$a,k))'^ = p(k) 
j^k 


GO 

ii) P(k+N) = Y_ ^(J.k+n)'^Q(j)i>(j,k-fn) 

j=k+n 

CO 

= Y *0.k)'^Q(i)«(j.k) = pW 

since Q(j+N)=Q(j) 

00 

iii) x'^P(h)x=:^ (^j,k)x)'^Q{j){^j,k)x) 

J-k 

>0 X 0 


(4.64) 


(4.65) 


(4.66.) 


(4.67) 


(4.68) 
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Slnce |#(j+k){ 0 if lA(k)| 0 V k *. P(k) is symmetric, periodic 

of period N, and positive definite. This completes the proof of the theorem. 

We note in passing that Theorem 6 is a special case of Theorem 8 when 

N = 1. 

{ 1 i i ) Linear Difference Equations with Variable Coefficients 
Given the difference equation 

= A(n) ^ (4.69) 

we can say very little about the stability of equation (4.69) for the general 
case of artitrary step varying matrices A(n). If the matrix A(n) can be 
represented as 

A(n) = + B{n) (4,?o) 

where A^(n) is either a constant or a periodic matrix, then in a number of 
cases we can develop sufficient, but not necessary conditions for stability. 

Theorem 9 

Given the linear difference equation 

Vl " ^0^^^ ^ ^ I\(n)I / 0 (4.71) 

where A^(n) is either a constant or a periodic matrix, 

% 

If i) V solutions of bounded as n “ 

f CO 

ii) I ||B{i)ll = b < « 
i=0 ° 


See Appendix 1. 
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Then V solutions of (4.71) are bounded for Kn . Before proving Theorem 9 we 
shall establish two important lemmas. 

Lemma 1 (Discrete Form of Bellman-Gronwall *s Lemma) 


n-1 

If 0(n) ^ C + i|j{i) e(i ) 


9(i)5»l'(i))C >0 J 

n-1 

Then 9(n) ^ C II [1 +ijj(i)] 
i=0 


(4.72) 


(4.73) 


Proof 

From (4.72) 


9(n) ijj(n) 

C+ [ Ijj(i} e(i) 
i=0 


< il'(n) 


(4.74) 


1 + 


6(n) i|j(n) 
n-1 

C + I i|^(i) e(i) 
i=0 


[1 + ii^(n)] 


(4.75) 


[C + r ijj(i) e(i)] ^ [C + I ^j;(i) 0{i)][l + i|;(n)] 

i =0 i=o 


( 4 . 76 ) 


n-1 n-2 

[C+ I i|^(i) 0(i)] < [C + I t|j{i) e(i)][l + ,j;(n-l)] 
i=0 1=0 


( 4 . 77 ) 


n n-2 

[C+ [ »jj(i) 9(i)] < [C + I i|;(i) 6(i)][l + i|i(n)][l+i;^(n-l)] (4,78) 
i=0 i=0 


n n 

Hence [C + J ij>(i ) 0(i) ] < C n [l+iji{i)] 
i=0 i=0 


( 4 . 79 ) 
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But 0n+i<C+ i ip(l)e(i) 
” ' 1=^0 


n~l 

6(n) < C n [1 + i|»(i)3 
1^=0 


(4.80) 
(4. 81') 


Lemma 2 


The product series S^: 


S„ = n (1 + V.) 

n i=0 ^ 


v^. i 0 


( 4 . 82 ) 


is convergent iff the series V 


n 

'n “ ^ ''l 
" 1=0 ’ 


is convergent. 


Proof: 


1) The product series (4.82) is convergent if the series L 


(4.83) 


= i An(l+v.) 

n 1=0 ^ 


(4.84) 


is convergent. This follows immediately from the fact that 




If Lim(Jln L ) = L, then Lim(S ) = e^ = S 


(4.85) 
(4. 86) 


2) We know that if is convergent, 0 as n-5“». let N be such that 

I 

for n >: N v„ ^ . 

n 2 
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Now 


? '^n = ''n(’ - T 


( 1 - 1 ) 


-) = - -y - lo •••) 


V? V? 


But 


An(l +v^) = (v^ - -J-+ •••) 

,2 

...) 


= VpO 


''ti ^ 

T'-'-T 


If 


then ^n{l +Vp) = Vp{l - ^ - ...) 


\\ < iin(l+vj < v„n+^+p+ ••■) 


<tn(H-v„) <|v„ 


Thus 


1) If I V. < 
i=0 ^ 


a) I V. < » and I v. < 
i=0 ^ k+1 ^ 


I And + v.) < p I V- < 
1=N+1 . ^ ^ 1=N+1 ^ 


(4.87) 

(4.88) 

(4.89) 


(4.90) 

(4.91) 

(4.92) 

(4.93) 

(4.94) 


( 4 . 95 ) 

(4.96) 


Returning now to the proof of Theorem 9, using equation (3.22), 


ii ) If I An(l + V- ) < “ 
i=0 ^ 

CO CO 

then I V. < 2 'I An(l +vd < “ 

i=0 ^ i=0 ^ 


(4.97) 
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Taking norms of both sides of equation (4.97) 

112^11 i l|X„|l l|c||+"i’ 11X^11 liX^’illllBCDIlIIXill 

But, by hypothesis i),||Xjj|| <M,, ||X|^|| ||xr|^|| < Mj |/1, 

llXnll Hell Y l|B(i)ll IlKjll 

Using Lemma 1 with C = M-j ||£|| , 0(i) = ||)i. j| and 
^(1) = Mg l!B(i)il 

we have 

ll2„IUM, Hell V (7 + M2l|8(i)ll) 

00 

But by hypothesis ii) I liB(i)l} = b < “ 

i=0 ° 

CO 

By Lemma 2, H (l + M, ||B(i)H } < d. < « 
i=0 ° 

il^il < M|dQl|c(l Vn 
Using hypothesis iii) we see that 

|IXj^|l<oo Yn 

Thus proving the theorem. 

Theorem 10 

Given the linear difference equation 


(4.98) 


(4.99) 


(4.100) 


( 4 . 101 ) 


(4. 102) 


( 4 . 103 ) 


2n+l “ ^<">4 


lA^Cn)) ^ 0 


(4. 104) 
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where A^(n) is either a constant or a periodic matrix. If 

i) A^(n) is a stability matrix, i.e,, solutions of = Aj^(n)x^ 
tend to zero as n ^ “ 

ii) llB(n)|l^bQ Jf n and sufficiently small. 

Then t^solutions of (4. 104) tend to zero as n ^ and the origin is Liapunov 
asymptotically stable. 


Proof. 


As before. 


n-1 


,-l 




( 4 . 105 ) 


Taking norms of both sides of equation (4. 105) 


n-1 


llx„IU ll^ll llcll + ,_M|X„|1 1 |x:Li|||b(i)||||.,|| 


i=0 


Using hypothesis i) [jX^H <. , 6 < 1 

Using hypothesis ii) equation (4. 106) becomes 


n-1 


||x„|| < M,||c||«'’ + M 2 b/J 6"-’-’ II X,. II 


i=0 


(4.106) 


(4.107) 


(4.108) 


Multiplying both sides of equation (4. 108)by 6 and setting 6(i ) =|| ll ; 
li£ll i using Lemma 1 


n-1 M„b 
< M,|(c|( n (1 

^ ‘ i=0 0 


(4.109) 


n-1 M„b MJ!c[l 

||x„||< M^llcll 6" .n^ (1 +4^) =-L— (S + Mgb/ 


5 


(4. 110) 
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Hence, if is sufficiently small, 


S + < 1 

z □ 


Ml Hell 

IlSnII i 


x„|| -^0 as n 
-fs " 


(4.111) 


(4.112) 


Therefore, the trivial solution of(4.,104)is L.A.S. 


Theorem 10a 

Given the same hypotheses as Theorem 10, we can prove the theorem using 
Liapunov's direct approach. 

Proof . Since A^{n) is a stability matrix, we know that there exists a sym- 
metric, positive definite, periodic matrix P(n) such that 


■ aJ(h) P(n+1) A^(n) " P(n) =“Q(n) (4.113) 

Q(n) = Q(n)^ = Q(rri-N) positive definite 

Let V = J P(n)x (4.114) 

n — n HI 

then P(n+l)x„^., (4.H5) 

Using equation (4.90) 

Vl “ Pth+l) V"»2n 

+ P(n+l)B(n)) + B(n)''' P(n+1) B(n)x^ (4.U6) 

^''n = Vl ■''n “ 


(4.117) 
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where S(n) = B^(n) P(n•^l) A^(n) + A^(n) P(n+1) B(n) 
+ B^(n) P(n+1) B{n) = S^{n) 


(4.118) 


Using (4. 112), equation (4. 118) becomes 


^ (4.119) 

1 

Since Q(n) is positive definite forvn, it is clear that by making j|B{n)f] , 
and hence the elements of B(n), sufficiently small, aV^ can be made negative 
definite. 

Hence for llB(n)|j sufficiently small, 

AV^ < 0 (4.120) 


V < V < V , • • 
n+1 n n-1 


< < ''o 


(4.121) 


Since V is positive definite and vanishes only at the origin, therefore V„ 
n n 

and hence tends to zero as n-s-«, and since is bounded above by V^, 

ll^ll is bounded for all n. Therefore the trivial solution of(4.l04)is 
Liapunov stable. 


Note: Let X(n) be the smallest eigenvalue of Q{n) and u{n) be the largest 

eigenvalue of S{n) in absolute value. 

Let r = Min X{n), r is positive, since Q(n) is positive definite, 
n 

s = Max p(n), we note that since S(n) tends to zero as b tends to zero, 
n ° 

s may be made arbitrarily small by making sufficiently small. Now 

4V„ < 

and hence by making b^ sufficiently small AV^ can be made negative definite. 
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b) Stability of Nonlinear Difference Equations 

(1) Stability of Explicit Nonlinear Difference Equations 

Theorem 11 (Liapunov-Poincare) 

Given the nonlinear difference equation 




(4.122) 




where A(n) Is either a constant matrix or a periodic matrix, 

If 1) A(n) Is a stability maL. ix 

„) 11„ UI(MlU=o Yn \ 

Hill -0 ||x!| 

111) ||£|| is sufficiently small 

then V solution of equations (4. 122) are Liapunov asymptoti cally stable. 

Proof . If A(n) is a stability matrix, then by Theorem 8 there exists a 
symmetric, positive definite, periodic matrix P(n) such that 


(4.123) 


A(n) P(n+1) A(n) - P(n) = - Q(n) 

Q(n) = Q"^(n) = Q(n+N) positive definite 


(4.124) 


Let P(n)x„ 

Vr 4i 


(4.125) 

(4.126) 


Making use of equations (4.122), 


Vl "" P(n+1) A(n})^ + x^(A'(n) P(n+1) f(x^^n)) 

+ f(^,n)"^{P(n+l) A(n))j^ + P(n+1) f(Xj^,n) 




‘-n- 


'-n 




(4.127) 
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Thus AV^ = V 

n n+l n 

= xJ(A^{n) P(n+1) A(n) - P(n))x^ 
iiJ(A^(n) P(n+7) f(x^>n) + f”^(x^,n)(P{n+l ) A(n))Xj^ 

+ P(n+T) f(^,n) 

— n — ' 'HI 

Using (4.124) equation (4.128) becomes 

= -xj Q(n)x„ 

xJ(A^{n) P(n+1) f{x„,n) + ,n)(P(n+l ) A(n)))^ 

+ P(n+1) f(^,n) 

Using hypothesis 11), lit(x,n)|j 0{|l2ii*^) as ||xlj -»• 0 
Hence a) ^ Q(n)xJ 'x- 0( Hx||^) 

b) iiI(A^(n) P(n+1) f(x„»n) ■*- f^C>^>n) (P(n+1 ) A(r\))x 0( ]|x 

c) P(n+1) f(>^>n) ^ 0(|j^|(^) 


as Il^ll ^ 0 

Thus for ilxj^il sufficiently small, the sign o'P is that of the 
.■ AV^ is negative definite. Hence, 


V XT < V < V , < 
n+l n n-1 




Thus if [|£[j Is sufficiently small 


aVj^ < 0 , 1/n 


and since Is positive definite and vanishes only at the origin 
fore 0, and hence llx^^II "^0 as n -j- 


(4.128) 


(4.129) 


iP) 


(4.130) 
first term 


(4. 131) 


(4. 132) 
there- 


Thus equation (4. 122) is 



t 


■ I 

■ -i 
■■1 
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Llapunov asymptotically stable at the origin. 
Theorem 12 

Given the nonlinear difference equation 


= [A^(n)+B(n)]x^ + 



where A^(n) is either a constant matrix or a periodic matrix* 

If i) Ap(n) is a stability matrix 
ii) ||B(n)|| is sufficiently small 


i i 1 ) 



llf{x,n)l 


= 0 


V^n 




(4.133) 


(4. 134) 


iv) Ii cll is sufficiently small J 

Then V solutions of equation (4.133) are Liapunov asymptotically stable. 

Proof . The proof follows along exactly the same lines as Theorem 10a and 
Theorem 11, and will not be repeated here. 


ii ) Stability of Implicit Nonlinear Difference Equations 


Theorem 13 


Given the implicit nonlinear difference equations 
Vl =A(n)x„ + f(x„,x„+,.n) 

[A(n)| i 0 

Vn 

||A(n)|| < “ j 


-0 — 




(4.13.5) 


where A(n) is either a constant matrix or a periodic matrix. If 
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•f) A(n) is a stability matrix 

i i ) 1 1 tn 






■’’’’i) ll^il sufficiently small 
Then V solutions of equation (4. 135) are Liapunov asymptotically stable. 


{4.136) 




Proof . Since A(n) is a stability matrix, then by Theorem 8 there exists 
a symmetric, positive definite, periodic matrix P(n) such that: 


i) ACn)"^ P(n+1) A(n) - P(n) ^ -Q(n) 
ii) Q(n) = Q^(n) = Q{N+n) positive definite 


(4.137) 


Let V = ^ P(n)x > 0 
n -n -n 


><n ^ 0 


“ 4+1 


(4.138) 

(4. 139) 


Making use of (4. 135), 

P(n+1) A(n))>^ + xJ(A'''(n) P(n'M)) .n) 

I^(i^>iinH.-j.n)(P(n+l) A(n))2^ 

+ (4.140) 


= >^(A^(n) P(n+D) f(x,^,Xj^^.^ ,ll)+x^(A'^(n)P(n+l)A(n)-P(n})^ 

A(n))^ (4.141) 
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Making use of equation (4, 137) 

Q(n)x„ + xJ{A^(n)P(n+l)) I(>^.x^+T,n) 

■*' i P{n+1) 3 ^) (4.142) 


From equation (4.135) 


Ill5n+ill + ll^ll <I|1 + Mn)|| llx^ll + ||i(^,x„+l.ti: 

lIVlII ■*■ ll^ll - lll(2n’^+l’'’)|l 


From (4. 136ii) 


where 


llI(i<n»-^^+r")|l ^ M2{6)(||)^|j + !l4+ill)‘ 

lli^ll + IliVi+lll ^ 


Mg (6) 'V'O(l) as 6"^0 


From (4. 144) and (4. 145) 

Millinll ,, ,, 

IliJn+lll + ll-<f|ll i 1 - ^ ^ 

i”*" lli^li < ff" is sufficiently small. 


(4.143) 
'(4. 144) 


(4. 145) 


(4. 146) 


Thus, if li^ll is sufficiently small, the first term in (4.142) is of 
order ||x^|| , while the second and third terms are of order ||j^|| , and 
the fourth term is of order |jXj^||^. Hence, if l|x^|| is sufficiently small. 
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the sign of AV^ is that of the first teriiis therefore 

AVj^ < 0 (4.147) 

Thus, if |i£|f> is sufficiently small, aV^ is negative for all n , there- 
fore 


" ''n ''n-l < • • • < V, < V„ (4, 148) 

Since is positive definite and vanishes only at x. = 0 , therefore 
-j* 0 as n ^ and equation (4. 135) is Liapunov asymptotically stable 
at X = 0, provided the initial data are sufficiently small. 


Theorem 14 (Liapunov-Poincare) 

Given the nonlinear difference equations 

(4.149) 

—0 — 

where A(n) is either a constant matrix or a periodic matrix, 

If i) there exists at least one unstable solution of the equation 
Vi = 

, !ll(2L,n)|i 

ii) lim = 0 

llxlhO llxll 

iii) |{£|| is sufficiently small 


Vn 


) (4.150) 


Then there exist unstable solutions of equation (■4..149). 
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Proof. Let 




«k = \ « 


-k 


Let Xk = 


(4.151) 


(4.152) 


Substituting into equation (4.149) 


= ®nll ®nll £'Vn>"> 


^0 


(4. 153) 


Now 


.-1 


n+l„-l 


e-|l A(n)e„ = R"^'X-;,A(r.)X„R 


1-n 


But = A(n)X|^ 


(4.154) 


e'L A{n)e„ = R 


(4.155) 


in+i = 


(4. 156) 




= b ; b-e‘’c ; a(y„.n) - 5“|, f(e V ,n) 


Since R = therefore from (4. 1501) r must have at least one 

eigenvalue of modulus greater than unity. 

Suppose that R is simple, and that the first k eigenvalues have modu- 
lus greater than unity, suppose that the remaining (L-k) eigenvalues have 
modulus less than unity. Since R is simple, there exists a similarity matrix 
T such that 
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r^RT = A 


wh&re |x^*| >1 i.£ (l,k) 
[Aj.|<l je{k+l,L) 


Let 

Then 


y = T z 

^ HI 




^ = d 

—0 — 


where h_ = T“^£{T^,n) 


Let 



1 

I— 

II 

H3 



i 

1 


P = 

Ik 

0 


0 

'^L-k 


''n = ^n% 


It will be observed that is sign indefinite. 

U =7* P 7 

V+l 4l+l -^+1 
Substituting from (4.145), 

InH “ ^ P ^ ^ A z 

— n-r I — n HI — n m 

+ 2^ A* P ll(^.n) -5- ]V^(2^,n) P il(z^,n) 


(4. 157) 


(4,158) 


(4.159) 


(4.160) 


(4.161) 


(4. 162) 
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AV = V - V 
n n+1 n 


= ^[A*PA - 


+ Ch*C^,n)P A^ + ^A*P h{^,n)) 






(4. l63) 


Now 


A P A - P = 

Since 
|X,| > 1 


ixjl < 1 


(1^,1^ - 1 

0 

0 

0 - 

i € (l,k) 



G e(k+l,L) 


(4.164) 


(4. 165) 


A*PA - P is positive definite Hermitian. Using (4. 150) iii , the first term 

n 

in (4.163) is positive and of order 11^11 > the second term is of order 
while the fourth term is of order thus for sufficiently 

small ll^ll i the sign of is that of the first term and is positive. 


AV^ > 0 for 11^1 j sufficiently small, (4.166) 

Since V., is sign indefinite we can define a set 

il: V^>0; ilzjl <5 (4.167) 

Clearly, the origin is a boundary point of S2. In O , > 0, AV^ > 0, 

therefore starting in n, z cannot approach the origin. Since V >0, 

0 







^ only exit S? through the boundary ||4,j| = 6; thus the system is unstable 


Theorem 15 




Given the nonlinear implicit equation 


Vl ' Mn)x„ + 


^ = c 

-iQ — 


where A{n) is either a constant or a periodic matrix. 

If i) there exists at least one unstable solution of the equation 


(4.168) 


4+1 


ii) lim 


f(x,y.,n) 


!|x|U|y[ho lixjl + [|y| 


= 0 


iT^n 


(4.169) 


ll^il is sufficiently small 

then there exist unstable solutions to equation (4,168) 

Proof . The proof follows along the same lines as that of Theorems 13 and 
14 and will not be repeated here. 

Theorem 16 (Liapunov-Poincar^) 

Given the nonlinear difference equation 

i<n+i = A(n))t + ffe.n) 


4 i 


(4. 170) 


-45- 


where A(n) is either a constant matrix or a periodic matrix. 

If i) the principal matrix of the linear difference equation 

4+1 = 

has an eigenvalue +1, or a pair of complex conjugate eigenvalues 
of modulus unity 


n) lim = 0 

llxll+o lixil 


Vn 


iii} ll^ll sufficiently small 


J 


{ 4 . 171 ) 


then the stability of equation (4.156) cannot be decided from the stability 
of the linearized equation. 


Proof . If we repeat the proof of Theorem 14, we see that in this case, 
(A*PA-P) is only positive semidefinite, having a zero eigenvalue corres- 
ponding to X = ±1 , or a pair of zero eigenvalues corresponding to |x[ - 1. 

Since' the matrix (A*PA- P) is only positive semidefinite, we see that 
the sign of depends on the terms in il(Zpin), Thus the stability is 
not determined by the stability of the linearized equations. 

Theorem 17 Theorem 16 is easily generalized to the case of implicit non- 
linear difference equations. 

Theorems 16 and 17 cover what are known as the "critical cases," that 
is, those cases in which the stability is not determined by the stability 
of the linearized equations. 
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5. DIFFERENTIAL EQUATIONS AND DIFFERENCE EQUATIONS 


(a) Numerical Solution of Ordinary Differential Equations 

As pointed out in the introduction, one of the more important sources 
of difference equations occurs in the numerical solution of ordinary differen- 
tial equations. 

Given the system of differential equations 
d_x 

= A(t)^ + I 

( (5.1) 


x(T:q) = c < t < T^< J 


we wish to approximate the solution of equation (5.1) by the solution of the 
di f fe ren ce eq u a ti on 






(5.2) 


such that y^^ ~ ^n+1 “ V n = 0,1 ,2 , • • • 

The natural requirements for the approximating difference equations are 
that for any function class of sufficiently differentiable func- 

ti ons 

1) They have a unique solution, 

2) This solution, at least for sufficiently small At^, should be close 
to the exact solution of equation (5.1), 

3} This solution should be effectively computable. 

These three points are examined in detail in books on numerical analysis and 
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will not be pursued at length in this note. 


(b) Numerical Solution of Linear Ordinary Differential Equations 
Consider the system of differential equations 

dx 

— = Ax + f(t) 


0 < t < T < 

' '■ ' r\ 


(5.3) 


X(0) = 
One technique 


c A = a constant matrix J 

for solving (5.3) is the use of the trapezoidal algorithm 


=^+1 ~ ^ 2 ^ '^n+1 ^ ^ 2 ^ -^n+1 ^ 

y^ = c At = VM 

Equation (5.4) may be written in explicit form; 

where 

= [I - fl [I + A |5-] 

■® = [I - A ^ 

Accuracy 

Let x^^-j = x(nAt) 

5e the local truncation error defined by 

+ 4) Vi 


> 


> 



(5.4) 


(5.5) 


(5.6) 


(5.7) 


Let e„ = (x„ - y„) be the solution error 
—n —n — n 


(5.8) 


Then subtracting (5.5) from (5.7) 


e . T = Ae + T , T At 
— n+1 -n — n+1 


(5.9) 


Now = x„ = 0 
— 0 — 0 ^ 


"1 


% 




= (A^-j + Tg) At 


e 

-n 




i.) 

—n 




J 


(5.10) 


If the matrix A is simple, there exists a similarity matrix T such that 


T"^AT = A 


and hence 


A = TAT 


-1 


(5.11) 


Hence If the homogeneous solutions of (5.3) are stable we know from, theory 
that X(A) must either be pure imaginary or have negative real parts. 


A= [I-A^f’[I + A|t] = t[I-A^]"''CI+A^] r'' 


(5.12) 


A = T 0 T' 


-1 


where 0 = [6.] 

J 


e- = (1 - X. f )-i 0 Xj f ) 


(5.13) 


If X. is pure imaginary, say X. = iw. , then 
J J J 
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If A. is complex with negative real part, say 


than 


ISjl - 


At,?- 


(“i T-) 


ni/2 


1 + tOjSjAt + (tOj 


< 1 


- + i(o- 

J^J 3 



(5,15) 


Hence the horTioganeous solutions of (5.6) are Liapunov stable by Theorem 5. 

In either case we have 

= T[0""V“'’t^ + e"'^T"’'T2 + T“’’t^] At (5.16) 

Let T"^^. = b. 

IlinI! ^ ilTlI [||Q""^b^ll + lib^ll ] At (5.17) 

But il0"-\ll = I |0,r-‘^|b^I< I ib^|= ll^ll (5.18) 

“K j=l J J j=i 3 -K 



(5.19) 


If T = Max|KI| 

k ^ 

then n^ir< ||Tll lir'^llnAt r 
but nAt = 5 T 

IKII iT||T|l ]|T-’||t 
< 


(5.20) 

(5.21) 

(5.22) 

(5.23) 


Now 
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^+1 


^n+1 


X- + 7j^ At + ^ + 

-n dt .*2 2 


dt^ 


df 


^ dt ^ ..2 2 ^ ..3 6 


W.3 

T' 

3 


dr 


dt' 


Substituting into {5.7} and using (5.3) 


^+1 


Atn-i r 1 n3.. .1 yi2 


dfn 1 d^f 5 
m , I — TIt ,J^2 




dt 


df 


d f. 


Then, provided \\>^\\ , ||f^H , ||^H , ][— “H are bounded. 


'-n’ 


dt‘ 


1:^+1 If - 


T = K2At 


as At ■> 0 


l!^illK^I< 2 Ar as At 0 

The trapezoidal scheme is second order accurate as At ^ 0. 
Application 

Consider the conservative dynamical system: 

Mx + Kx = 0 


^(0) = ^^x(0) = b 
where M - is positive definite 

j 

K = K is positive definite 
'x 


If 


equation (5.28) may be written 


(5.24) 

(5.25) 

(5.26) 

(5.27) 


(5.28) 
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2(0) = c = (b) 


where A - 


0 I 


K 0 


If A is simple, there exists a similarity matrix such that 


r .. 


T“^AT = 


lOJn 


- 10 ), 


10)^ 


-Vlir 


The trapezoidal difference equation corresponding to (5.29) is 


^.,1 = Wp T ^^Vl 




4 




w„ = c 

— o — 


Alternatively, 


^+1 ^ ^ ^ 2 ^^+1 




^+1 “ ^ 2 ^ ^^=^+1 ^ 


J 


From {5.33} we see that 


j ij+i H y„+i + S' 4i Vi 


(5.29) 


(5.30) 


(5.31) 


(5.32) 


(5.33) 


= constant 


(5.34) 
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Thus the difference equations (5.32) or (5.33) conserve energy in exactly 
the same way as equation (5.28) whose first integral is 

^ 1 )Jmx = constant (5.35) 

From equation (5.32) 


% = c 

where [I - ^ A]""* [I + ^A] 

Using (5.31) 

T[I - ^ A]“^ [I + ^ A] T"'*' 


= T0T“^ ; 


0=[e^-] 




Vj 0,2N) 


Hence |Xj(^)| = 1 


(5.36) 

(5,3;) 


(5.38) 


(5.39) 


(5.40) 


Thus the eigenvalues of u?all have modulus equal to unity and equa- 
tion (5.36) is Liapunov stable. This property is exemplified in t' fact 
that the energy is conserved. 

Using equation (3.5), the solution of equation (5.36) is: 





where u?" = T T“^ 


(5.41) 

(5.42) 
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1 + 1 4^ tCj 

, -T- a ■ = ® 

■> 1 - 1 T “j 

ji) jAt 

, -1 

where A^. = tan * . 2 


Let R . ^ 


^ = T e 


^-1 


(5.43) 


(5.44) 


(5.45) 


The solution of equation (5,29) at t^ = nAt is: 


— n 


^ V T T-' 


(5.46) 


We see that the solution of the differential equation (5.29) and the cor- 
responding difference equation (5.36) have the same structure, however, in 
general the time dependence is different. 

Period Error 

Let T? = ^ be the period of the j'*'*’* mode of the difference equa- 

j 

tion. Let T. = 2ir/to. be the period of the mode of the differential 

J J 

equation. Then 

T^ - T.. 


is the period error of the j'th mode of the difference equation. 
We note that 


(5.47) 
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L1m S3. = Lim 
At^O ^ At^O 


1 

At 


’-<“j 


^'i2 
2 ^ 




(5.48) 


Thus in the limit as At - j- 0, the period error vanishes and equations (5.45) 
and (5.46) are identical. 

From a practical computing standpoint, we cannot let At go to zero. 
While {(jj.-At) can be made acceptably small for the lower modes of a complex 

•J 

structure, it is not possible to make (w.At) small for the highest modes. Thus 

ij 

by making At sufficiently small, equation (5.45) will give an accurate repre- 
sentation of the low mode behavior, however, higher mode behavior will not 
be accurately modeled. In most problems in structural dynamics, only low 
mode behavior is of real significance, therefore if high mode behavior can 
somehow be suppressed, equation (5.45) will give a reasonably accurate repre- 
sentation of the response of a complex structure. 


Methods Proposed for Suppressing the Higher Modes 


(i) Use of Viscous Damping 

By analogy with continuous time systems it might appear that the use 
damping could be used to suppress the higher modes. As we shall now show, 
the method is ineffectual in suppressing the higher modes in discrete systems. 
If in equation (5.28) we add viscous damping, then the equation becomes 


Mx+Cx + Kx = 0 ^ 


)^(0) = ^ x(0) = ^ 

M = > 0 ; C = >: 0 ; K = > 0 

If Z = equation (5.49) may be rewritten 




(5.49) 
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A z 
dt 


z(Q) = c = I 


(5.50) 


where A = 


-M"^K "M"^C 


If A is simple, there exists a similarity matrix T such that 


T"^AT = A = 


-co^^l+iw 






-W]5i-i'WnVl-Ci 


-WnSn-'^n 




(5.51) 


The trapezoidal difference equation correspondi ng to (5.50) is: 


Vi 


w„ = c 
-o — 



= A w 
-n 


(5.52) 




where 


A = TQT 


-I 


0j= 


Hi 

UUj.Cjf - 1-o.jf l-?2 




(5.53) 




9. may be expressed as 


-56" 


0 * = p - s 
a ^3 






1 - f ff 

1 + + {Wj 


A(j)- = tan 

u 


-1 


3 3 

(o.At /l-S 

J— 3 


{5.54) 




oj.At p 
1 

2 


l-{4r-) 


For the lower modes p. 1. As mode order Increases, p. decreases initially 

J J 

then starts to increase again. For the higher modes p- tends to unity. Thus 

J 

we see that viscous damping is not effective in suppressing the higher modes 
( i i ) Use of Algorithmic Damping 

If equation (5.32) is modified to read 


0 < a < 1 


HI 




^ = c 
~o — 


Equation (5.36) now becomes 


=Jl w„ 
-n+1 a -n 


(5.55) 


where 


= [I - (l-a)CitA]“’’ [I +aAtA] 

UC 




o? = T 0 T 
a 


-1 


where 0 


a-i 


1 + At tfnij 
1 - At(l“a) (1). 


J 


(5.56) 


(5.57) 
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Hen ce 


9 . = p e 


lA^aj 




% 


a-i 


1 + (oAtw-)^ 

*J 

1 + (l"a)AttOj)^^ 




(5.58) 


A^ = tan 

Oi - 

j 


•1 


oj.At 

J 


1 - a(l-a)(w^At)' 

J 


We note that 


i) 


P. = 1 
J 

when a - 

1/2 

ii) 

0 < a < 1/2 

p. < 1 



Thus 

in case ii) p. =i 

1 for the 

lower modes. 

while Pf =i for the 

J 1-a 


higher modes. Unfortunately, if a 1/2, it is easily shown that 


ll£„ll < '=3lj-ffl|0(At) -H<,K2 O(At^) (5,59) 

Thus, If a i , the modified trapezoidal altorithm' (5.54) is only of first 

order accuracy as At 0. 

(i i i ) Use of Temporal Filtering 

1 1- T 

Let V = f 2iv + w , J 

— n 4 — n**l m — n-l 

where 

^ 21 


(5.60) 

(5.61) 


(5.62) 
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Now, for the trapezoidal algorithm 


Jl ^ T0T" 


where 0. = 


1 + it 04 


1 - ito. ^ 


(5.63) 


(5.64) 


1 T[0 + 21 + 0"'']T'’ 

1 + iw 1 - iw 

J_ T bi e_ + ? + 

' 1 - i«.j f 1 f 


i, = T 
— n 


3 7 "^ vf 

I + ^)2 J 


T-1 


(5.65) 


(5.66) 


(5.67) 


We note that Srr-o - 1 for the low modes and tends to zero for the 

high modes. Since ^ is second order accurate as At + 0, tJis filter, which 
is also second order accurate, still retains second order accuracy. Thus, 
unlike algorithmic damping, the use of the temporal filter does not affect 
the accuracy of the computational scheme. 

There exist many more sophisticated algorithms for solving problems 
such as equation (5.28), however, the author's experience has been that for 
linear problems, the trapezoidal algorithm with post-filtering does as good 
a Job as the more sophisticated schemes when applied to large complex struc- 
tures , 

(c) Numerical Solution of Nonlinear Ordinary Differential Equations 


Consider the system of nonlinear differential equations 


^ = A 21 + I(^) + £(^) 

x(0) = G 


0 < t < T < “ 


A = constant matrix 


(5.68) 
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= c 


One technique for solving (5.68) is the use of the trapezoidal algorithm 

^n-M = ill f - ^ Sn+1 

= c At = T/M J (5.69) 

Equation (5,69) may also be written as 

y„+l = + f(Vn) + Vi + Sn) 


(5.70) 


Accuracy 

= x(nAt) 

l- 2 t be the local truncation error defined by 

^n+r * i<4+l> i(2n> * Vl + 

(5.71) 

Let ^ he the solution error (5.72) 

Then subtracting (5.69) from (5.71) 


Vl " A ^ ^2^ Vl 


+ il(^) - Vl 

■^) bounded Vn e(ljM) 

ii) IIKx) ~ i(z)l! 1K|!2L-Z11 V 2 ^ hounded 

iii) f(>^) continuous and continuous first and second 
partials 




(5.73) 


> (5.74) 




-> 60 ^ 


Then 


llVlIi ^ <1 H|A||f]||e„|| + llAllf ||e„^lll 


2 IlSn+ill I' 2 IlSn^l ’’’ 


At 


Is 


%i+l" - 


Cl + (||A|1 +K) f lIln+lllAt 

" II + -n 


(1 - ( A +K)-^)'"^" ■ 1-(||A||+K)^ 


1^1 


= 0 


Thus [}e^ II £ 


IlSoll < 


IX] II 


l-( IIAII +K) 


At 


l + C llAll +K) 


At. 


Ll-( ||A||+K) f-J 


!t] I! IW 5 


At 


1 - ( llAll +K) ^ 


e 

i_ni 


< [f""’ 111 , 


t > 1 -""^ 1 |t ,11 t 


lh„ll ] — 
n 1 - { 


At 


+ K) 


At 


where 


1 t ( llAll + K) ^ 
1 - C llAll+K) ^ 


If 



At 

1 - ( ||A|| +K) 


But 


c];-l 


U( ilAll -HQ f- ^ ^ At( llAll -FK) 

1 “ ( ||A|| +K) ^ I - ( ilAli +K) ^ 


(5.75) 


(5.76) 


(5.77) 


(5.78) 


(5.79) 
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• 
■ ■ 

llg„ll |-7T — - 0 

^ ||A||+K 

(5.80) 

Now 

= 1 + y + ^ 0 < 0 < 1 

(5.81) 

• 

** 

(1 +y) i e'^ 

(5.82) 


1 - y >: 1 - 2y + i (2y)^ for y ^ 

(5.83) 



(5.84) 

. 

-- 1 - ^ f ^ 

(5.85) 


2 3 “ ® 

1 , 4. . z'^ „-0Z 

1 - z T e + e 


* 

.-■*- - -1 ^ n^y 

(5.86) 


7 — V — 2 ^ 

^ 1 - 2y + 2y 

* 

•* 

iix < Jy 

i-y -® 

(5.87) 

Hen ce , 

If (l|A|l + K) 1 



l + (||A|l+K)^ |(||A|| + K)At 

i> = £ e'^ 

1 = (IIAII +K) ^ 

(5.88) 

• 

»• 

|( llAll +K)nit 

^ ,||A||+K 

(5.89) 

But 

nAt = tj^ <. T 


• 

•« 

,, ,, , |( IIAII ^-K)T 

IlSnlU ,, „ e 

||A|1 +K 

(5.90) 


Returning to equation {5.71), 
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-^+1 " ^ dt 2 


-ti 


i!Vl) = f(Xn) + J(x„)#At-f... 


(5.91) 


da 

-n 


Vl " Sn+^At+etc. 


J 


Then using (5.71) and (5.68) it may be shown that 


llln+lll <K3(4t)' 


as At ^ 0 


iisnii i 


K,(4t)^ 


|( IIA1I+ K)T 


llAll + K 




as At 0 




(5.92) 


(5.93) 


Thus the trapezoidal scheme is second order accurate; unfortunately, unlike 
the situation for linear systems, the trapezoidal difference equations for 
nonlinear differential equations are not guaranteed to be globally stable. 


Application 

Consider the conservative dynamical system 
M }[ + K X-.+ f (>^) = 0 
^(0) = a ^(0) = ^ 

M = > 0, K = > 0 

X(x) = ^U(x.) U()^) >0 21 0 

Equation (5.94) has the first integral 


(5.94) 
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^ X + ^ X + U(x) = const. 


If 


( — j equation (5.94) may be written 

X 


dz 


A z -)-£{z) 


z(0) = c 


where A = 


0 I 


-M'^K 0 


a(z) = 


“M'Vcx) 






If A is simple, it has the representation 


A = TAX' 


-1 


A = 


1 UJ n 


-lOJ. 


iw- 


'IWo 


etc. 


The trapezoidal difference equation corresponding to (5.94) is 


w„+l = «n + ^ w„) ^ T a(Wn>> 


At 




zD. 


Wo = ^ 


J 


A1 ternati vely. 


" in ^ T- 'in+1 ^ in> 


in+l = in - T 


J 


(5.95) 

(5.96) 

(5.97) 

(5.98) 

(5.99) 

(5.100) 
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From (5.100) we see that 

+ 1 (Vi (f(y„+i) ■'■Kin)) = ® <5.101) 

which can be written 

f i^+l My„+, + } Ky„^, + = constant (5. 102) 


Where ^ 


(5.103a) 


We note that if are bounded, then as At “S- 0, {Z-j+i "^-j) ^ that 

(5.103) becomes 


n+1 

L1ni U , = Lim ^ 
At^O ’ At-^0 1=0 


J (Vi'+T - (KXi+i) ■' 


(5.103b) 


■^U(y^) = U(^^^) 


In this case, equations (5.102) and (5.95) are identical. In general, for 
At finite, cannot be guaranteed to be positive, in which case equations 
(5.99) and 5.100) are not guaranteed to be stable. It should be pointed 
out, that equation (5.99) can be rewritten as 


!!n+l = -SWn ■^^(£{w„+i) •^ al%)) 
% = £■ 

where .j! = [I - ^ A]"'' [I + ^ A] 

3= [I - f-A]'’ ^ 


( 5 . 104 ) 
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Si nee thf. eigenvalues of A are pure imaginaryj the eigenvalues of J! all 
have modulus unity, thus equation (5.104) is one of the “critical cases" in 
Liapunov stability theory as discussed in Theorem 16. 

To illustrate these problems let us consider the following scalar 
problem: 


X + f(x) = 0 

f(x) = X |x| ^ 




(5.105) 


f(x) = sgn X + y(x - sgn x) |x| > 1 


Equation (5.105) has the following first integral 
^ x^ + F(x) = const. 

F(x) = j lx| < 1 > (5.106) 

= 1+ |x| +f {|x(-l)2 ; |x| > 1 

Since F(x) >0 x 0, equation (5.105) is globally Liapunov stable 
with respect to the origin. 

The trapezoidal difference equation corresponding to equation (5.105) 


- y„ = -n(f(y„+i) + f(y„)) 

which may be written 


At/2 I 


(5.107) 


^(%i > 


+ y 


n+l 


= z = y 


+ 2nyj.^ “ n f(y^) 


(5.108) 
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Since f(y) is piece-wise linear, equation (5.108) may be inverted to 


01 ve in terms of 


Thus = g(z„) 




n+r 


y„-i-n y„ 3i\) 
n 


(5.109) 


J 


1 ? 

Where g(z) = — - „■ z for |zj 5 1 + q 

1 +n 

= sgn z + — p- (2 - (1+ri^) sgn 2 ) for |z| > 1 +n^ 

1 + un 


> (5.110) 


|g{z)| < }zl ^2 


Vl " '^n+l 


n+1 


'^n+1 


(5.111) 

(5.112) 


Vl = = "n-1 


(5.113) 


Equation (5.113) may be written in several different forms, two of which are 
given below: 

.2 


a) 


Vl - ^n-1 


4n 


1 +n 


2 


° l^nliO + n^) 

g 

■ = n +n^) sgn z^^ + hLII D- ) _ (i 

1 + pn 

r 

for |Zj^l > 1 + n' 


) (5.114) 


J 


From (5.114) we see that if y > 0, sgn(z^^.| - 2Zj^ + z^^.j) = - sgn z^ (5.115) 
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Thus the sign of the finite difference curvature is always opposite 
to that of the displacement, thus the solutions are always oscillatory. 


b) , - 211^ z . z„ , - - k,(z ) 

(i+pn^) ^ (1 +yn^) ^ ^ 


•<5(2) = — ^ 2 for Iz! £ 1 + n^ 

1 

" sgn z for I2I > 1 + 

Equations can also be written in the first order form. 


(5.116) 


^+1 + Ji(9n) 


/ Zn ^ 

0 = 

Jl = 

^2(1 -yn^) 
1 +yn^ 

"1 



1 

0 __ 



^ 4n'^(l-u) 
1 + 


f<2(Zn) 


\ 


/ 




The matrix^ has eigenvalues 


i 

1 + 

= e±1<l> 


cos (|) = ( 


1 2 
1 “ un ^ 

1 +yn 



(IjlHQ 

1 +yn 



(5.117) 


(5.118) 


Using equation (3.5) with y 7^ 0 
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0 = fin'll 

•^1+1 ^ 

IJ 

0 + 1 




But 

TAT"^ 

;. 1= 

where 

PfI = 1 

i == 1,2 


,n-i 




(5.119) 


(5.120) 


-1 


< IITII IIT-'II lle^ll + ||T|1 ||T-1|| J^||k(0.)H (5.121) 


But lli{0i)il5V 


iML 


1 +yn'^ 


11^+1 II < K( + W 


1 -u 


1 +yn‘^ 


n) 


(5.122) 

(5.123) 


Thusj even though (5.117) may be unstable, it is only weakly unstable, with 
at most linear divergence. Now using (3.5) with p = 0, in this case 
A.j - 1, i = 1,2 

1 1 
0 1 

n 


A = T 


r“l 


(5.124) 


llVlIlillUl llT’hl Kll (2+n) + |T|1 (2+0 


-1 


< K[ lle^ll (2+n) + 4n^(2(n+l) + (5.125) 

n 

Thus, even in this case, the 0(n ) as n as it is only weakly 

■^1 

unstable. 

Equation (5.117) defines a continuous mapping M(*) such that: 

Vi ’ 

Therefore, by the Brower Fixed-Point theorem, there exists at least one 
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fixed point, or equilibrium solution. In the case of equation (5.117), it 
is easily seen that the only solution of 

i* = is i* = 0 (5.127) 

1 / 

Next, let M , k an integer, denote the mapping H applied times. Then 
there may exist a sequence of distinct points ^"'^(1),^*(2) such 

that 

m = l,2,---,(k-l) 

0*(1) = M*^(G*(D) 

Clearly, this sequence constitutes a periodic solution of period k. 

Stability of Periodic Solutions 


Let 0*{in) • tn ^ (l,k-l) be a k periodic solution of equation 

if 


(5.117). 


Let ^(m) = £*(m) + 5^(m) 

(5.129) 

Then 5£(nW-l) = M Q(_0*(rn)) 6£(m) 

(5.130) 


Provided £(m) and ^*(m) are on the same piecewise linear branches of ]<(£)• 
Thus, for 5^(k) small, but not infinitesimal, and ^*(m) not on a corner of 

m 

<S0(mM) = [ M^Q(0*{i))]5e{l) (5.131) 

K 

Hence 60(K+1) = [ I M fi(e*(i))] 60(1) 

1=1 

= 50 ( 1 ) 


(5.128) 


(5.132) 
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If U^-(A[^)| < 1, i€(U2), the piriodlc solution is stable. 

If U^-(A|^)| > 1, i C{1»2), the periodic solution is unstable. 

If lX^.(A|,)| = 1, iC'djS), the periodic solution is stable, provided Aj, is 

simple, Othewise unstable. This is a property which is special to piecewise 
linear systems. 

Example 


If in equation (5.116) we set y = 0 and n = 1, we have 

z.t-2z +z t= ~4ko(z ) 

11+1 n n-1 '^2^ 

k2{z) = |- for jzi £ 2 } 

= sgn 2 for |zl > 2 

(a) If ]2 q| ,|Zi 1 < 2 , then 

Vl ^n-1 = ° ^n^l = -Vl 

"4 " ^0 ’ ^5 = "l 

t 

^2n+l = 

Since |2^| ijz^l < 2, it follovus that |zj <2 Vk 
From equation (5.117) with y = 0, n = 1, 


^+1 


0 -1 




(5.133) 


(5.134) 


1 0 


(5.135) 
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L- 


0 -1 

1 0 


m^q(0 (D) = 
of {< 2 ( 2 ). Thus 

6l(K+l) = Aj, 6S(1) 


Vi , provided 0*{i) is not close to a 




\ = 


.n 


0 -1 

1 0 


K 






Thus as long as the initial perturbations are small, not necessarily 
finitesimally small, the periodic solutions are stable 


(b) If 12 q|,|zi| >2, there exist periodic solutions with |Zj^I > 
In particular. If = N+1 , z-j = 3N-1 , then Tj,^ = 2{N+1) 

Proof . If Z|^ > 2 , equation (5.133) becomes 

Vi - ^ \--\ = 

With Zq = N+1 and z^ = 3N-1 

= (i+N) + 2n(N-n) 

Thus Zj^_-j = .3N-1 = 2-j ; Zjyj = N+1 = z^ > 2 

Vi ' -(N+n = -Zp < 2 

If z^< -.2, equation (5.133) becomes 


corner 


(5.136) 


in- 


2 /n. 


(5.137) 


(5.138) 


(5.139) 


z.t-2z +z 1 
n+1 n n-1 


= +4 


(5.140) 
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= -^1 

Thus for N < n < 2N+1, the solution for 0 < n < N is repeated with the nega- 
tive sign. 

•^(N+2) " (5.141) 

Thus there exists a periodic solution with initial data = N+1, 2 -j = 3N-1, 

with [ 2^1 > 2 Vn and period Tj^j = 2(N+1). Clearly, there exists an infin- 
ity of such solutions. 


Stability 


Since each point ^*(k) satisfies the condition ]z^| >2, each 
iH (9*(i)) is the same. 


Now 




2 
1 

= T 


2 

1 

<• 

-1 

0 


-1 

0 


= T 


T 


-1 


T' 


-1 


(5.142) 

(5.143) 

(5.144) 


Thus, though |A-(A|^)| = 1 i/ i 6'(1,2), Aj^ is not simple, hence the periodic solu- 
tions are weakly unstable, and will grow until ^(n) = ^*(n) + (S9(n) reaches 
a corner of k 2 (z), at which point the nature of the stability will change. As 
shown in (5.125), the global rate of growth is limited to 0(n ) as n 
which is still a rather weak type of instability. 


Globally Unstable Solutions 


Hughes ( i ) and others have exhibited numerically the 
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instability of equation (5.107) and hence of equation (5.133). 
bility of equation (5.133) can also be exhibited analytically, 
equation (5.133) with prescribed initial data 

Vl ■ ^ ^n-1 = 

k 2 (z) = z/2 for U1 1 2 

= sgn z for \z\ > 2 


= 1.5N + 2.5 


^ 3.5N + 1.5 


Since z„,Zt > 2 
0 I 


Zj, = A + Bk - 2k'^, provided z^^_^ > 2 
Using the given initial data 


- 1 ,2, 3, etc. 


= A ^ 1.5N + 2.5 


z^ = A + B-2 = 3.5N + 1.5 


t 

B = 2N+1 J 
Zj^ == (1.5N + 2,5) + (2N+1 - 2k)k for k <. N+2 


= 0.5N + 1.5 , Zj^^2 " -0-5N + 3.5) 

To determine we use equation (5.145) 


^N+3 " + 1.5) + 4 


= -(3.5N + 4.5) 


Since, ^n+ 3 minus two, we can write 

^N+2+k ^ "^^1 ^ ^ ^ 

where A-j = "^|\j 4.2 0.5N + 3.5) 

^N+3 " "*^^1 B " 2] •• = 2N+3 

Zfj+2+k " ^ * k<n+2 

^2wH-5= 0-5NH-5.5) = (z^t3) 


The insta- 
Consider 

V 

> (5.145) 

i 

(5.146) 

(5.147) 

(5.148) 

(5.149) 

(5.150) . 

(5.151) 

(5.152) 
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To calculate Z2N+6 '"eturn to equation (5.145), from which 


■2N+6 


(3.5N + 8.5) = -t- 7 


Thus at the end of a complete cycle 
T^ = (2N+5) 


and = Zq + 3 

^T,+l “^1 


are the initial data for the 
next cycle 




At the end of the next complete cycle 


Tg = 4N + 14 


and Zt = z„ t 2 X 3 


At the end of the complete cycle 
= k(2N + 3 + 2k) 




and Zt .= z + 3k 

Tr “ 


Zy = Zn + 7k 
‘k+1 ^ 


are the initial data for the next 
cycle 


J 


Returning for a moment to equation (5.101) with n = 1 > then 


I -i-n) (yn+l = ° 

If "^n+l ’-^n greater than 

then (5.157) becomes 


(5.153) 

(5.154) 

(5.155) 


(5.156) 


(5.157) 

unity. 
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r ‘^n+l -^n> -yn>'=S" Vl =9" ^n’ “ ° 


5-<^n+l-^> 'IVT'-Iyni= “ 


(5.159) 


From equation (5.106) with and of the same sign and both greater 
than unity, then, 


2 ^^n+1 ■*■ *^^^0+1^ ** ° 


(5.160) 


But F(x^^^) - F(x^) = - jxj 


2 ^^n+1 “^n^ Vl l"l^ni " ° 


(5.161) 


Thus the trapezoidal algorithm preserves the energy identity (5.161) if 
y^^^ and y^^ are both on the same nonlinear saturated branch. If y^^^^ and 
y^^ are both on the linear branch of the curve, energy is again conserved. 

If and y^ are not on the linear branch or not on the same saturated 
nonlinear branch, then in general energy is not conserved. 

Returning to equations (5.109) and (5,112), 

W‘ 9("n' 

Vl " - ? "n-1 - 

Since IZj^l > .2 ly^+i 1 > 1 > thus we need only look at the "energy” at the 
beginning of each half-cycle to see how it is growing. 

" F^k+1 l^k+ll 

(5.163) 

IVli = l"n - "ni = l"nl ” ^ 


(5.162) 
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From (5.156) 

+ [z^ + 7kj - 1 (5.164) 

As k 

O(k^) (5.165) 

But from (5.155) 

= k{2N+3 x2k) . 

As k ^ , T|^ O(k^) (5.1B6) 

Therefore, combining (5.165) and (5.166) 

Ej^ ^ Tj^ as k « (5.167) 

Thus confirming analytically what Hughes and others had obtained numerically. 

Equation (5.116) was carefully examined for y = 0 and n arbitrary; 
nothing essentially new was learned, except that even for n very small, but 
not zero, weak instability will still occur if the initial data are large 
enough. 

Equation (5.116) was carefully examined for jp] > 0, n both arbitrary,' 
for Ip| sufficiently small the system behaves in very much the same way as 
for p = 0. It is true that the system appears to have bounded solutions, 
however, the bound is of the order of (l/p) and hence can become very large 
for certain ranges of n. 

Since the results of this section are essentially negative, we shall 
not report all the work that was done to investigate the effect of nonzero 
p, the effect of small n and the effect of damping. Instead, we refer the 
interested reader to the PhD thesis of my student, B. 0. Westermo (2). 
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Algorithms which Conserve Energy 


As shown in the last section the trapezoidal algorithm, which was 
found to he very useful for linear problems, can for a- certain class of non- 
linearities lead to weak instability. In this section we shall look at several 
new algorithms which conserve energy. 

Consider the conservation nonlinear differential equation 


d X 


dt 


+ f(x) = 0 


0 < t < T 


where 


and 


32 


^ " I > 0 X 5^ 0 


xf(x) >0 r 0 




(5. 168) 


The system (5, 168) has the first integral 


-f F(3 c) = const (5* 169 ) 

and since F(3c) is positive definite, {5. 168) is Liapunov stable with respect 
to the origin. 

If we rewrite equation (5, 168) in the form 


d X ^ dF n 
dt 


This immediately suggests the algorithm 

At /a , • \ 

s'Cy^+i)-F(yJ 


y , - y = - At 
^ n+1 ^ n 


^^n+1 ^n^ 


(5. 170) 


(5. 171) 
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Cross multiplxcation immediately shows that; 




Or on summing 




= const 


If in (5. 168) 


f(x:) = (1) X + g(x) 






{ju + g(x) > V X ^ 0 


J 


( 


g(x) continuous with continuous first and second derivatives then (f 
becomes 


At 


yn+l - = T 

yn+l'^n = - ^ ^(yn+I+^n) ' 

which may also. be written in the form 






At 


5^n+i =:si. +-r +xten+i’2:n) 


(: 


5. 172) 


» 173) 


5. 174) 


.171) 


. 175) 


. 176) 


where 
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A = 


u)^ 0 




n 


-21J 


%'■ 

Ac ni racy 


/ 


At 


G(yn+i)-G(y^) 


\ yn+l-yn 


Now 






^ g(yn)+g \7j — 2 — g (5) ^ 






Ty 


n+l ^n' 


§ - 9y^+(l-e)y. 


n+l 


0 < 0 < I 


Thus if y = y are hounded Vn 
^n’ -^n+1 




< M Vn 


yn+l -yn 

« . 

Similarly if hounded Vn 

11 ^ ^2!:n+i'^^All ^ 

Hence from equation (5. 176) 


n+l~— n^^ ^ (N+M)At 0(At) for At small 


As before, we define the truncation error T 

—n+l 


(5.177) 


^5.178) 

(5.179) 

(5. 180) 

(5*181) 

(5.183) 


, by the equations 
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&.+r ^ +X(5^+i.£n) +ln+l^* 




(5.. 183) 


Using equations (5. 174), (5. 178) and {3. 183) it is readily shown that if 


are bounded for all n: 


lT,.^lll'-0(An 


(5.184) 


Defining tlie solution error (5.176) is subtracted from 

(5. 183) 


At 

^0+1 =£n+T^(ln+l+®n>+X(5„+i.^)-X(S;n+l-Z„)+I„+l^‘ (5.185) 

Il2„+lll lllj llA|l(l|fjn.i(l+iyD + lb<.(5n+l>£a)-X(5:n+l'2^>'' 


+ (S-186) 


Using (5.178)and (5.182) 


provided thai z , are bounded Vn where 

— n — n+i — n — "n+i 

^n+l ~ 


(5.187) 


|fi - Ma3£(ft^^^-t-l|Tj^_^jll) ~ 0(Ar) as At 0 


(5.188) 


Using (5.186), (5. 187) and (5.188) 

A +(W +K) “If \ 


(5.189) 
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If ^ =0 tlien as previously shown 


1I.J1 . , 


exi 


Wl +K 


Thus for T fixed 


lie 1 s K,j6 ->^0(At^) as At -♦ 0 


n 


(5*190) 


(5.191) 


Thus tills algorithm is second order accurate; in addition it is Liapunov 
stable witli respect to the origin. 

We note in passing that equation (5, 178) may also”be Vi?x':t-+-=*n 


G(y„ ,i)-G(y ) , 

■ — = 2 + 
n+1 n 




IT 


^^n+1 


5. = y 6. -f(l “9.)y 
1 ^n+l 1 ^ 


0 < 0 . < 1 
1 


(5.192) 


Thus if y^, ^n+1 bounded Vn, then for At tending to zero. 

G(y ^i)“G(y ) . , 

- Z [g(yn4-i)+e(yJ3 + M(At)^ (5.193) 

Thus this present algorithm may be considered to be a modified "trapezoidal" 
algorithm. 

Effect of Viscous Damping 

If viscous damping is added to equation (5.168), tlien using (5. 174) 
we have 

X + 2zx + (JU^x + g(x) =: 0 , z > 0 (5 . 194) 

Let 
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y(x,x) = + 2G{x)] 

V(x,x) > 0 if X, 3 £ 0 


(5.195) 


V(x,3c) = XX + ax^ + axx ^ 


«J XX + 2z XX + g(x)x 


Using equation (S. 194) 

V = -a(x^ +U)^x^ +xg(x)) 


But, from equation (5. 174) 

2 2 

UJ X + xg(x) >0 Vx 0 
V< 0 


(5.196) 


(5.197) 


(5.198) 


Hence V is a Liapunov function and the system (5. 194) is Liapunov asymp- 
totically stable at tlie origin. 

Consider now tire discrete form of equation (5.194) 


Let 


Wyn=-2-(y„+l+yJ 

i+r^n " ■■^[^®(yn+l+yJ'""^(yn+l+^nW - 

^n+1 


■^11 = 1 * 2G(y^)] 

is positive definite and vanishes only when y = y = 0 


V 


n+l =2C(yn+l+^yn+l)^ + 


n+1 


m+l- 


(5.199) 

(5.200) 


(5.201) 


(5.202) 
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AV = V -V 
n n+1 n 


+ + G(y^^j)-C3{yJ 

Using equations (5. 200) 




“ At 


(5.203) 


^^^n+1^ ^^^n^ zAt ,• • X At 2 

~y~r^ ^ ^W^^n^ - T (yn-f-l+^n^ 


X 


[■M (yn+i-yn) + ^(yn+i+y^)] 


+ i [z^-KU^][y^^j^ - y^] + - alY^) 

• 2 


(5.204) 


AV = -zAt 
n 




yn+I-ya 


(! 


If G(y) is an even monotone increasing function say 


G(y) = G*(y^) 

❖ 2 

G (y ) is also a monotone increasing function. Hence 


(5.206) 


(yn+l>-‘^ (yn) ^ . 
2 2 ^ ° 


(5.207) 


n+1 


AV„ = -zAt 
n 


/ yn+(+y„ \ , ''n+u ’ ^^n 

l— +(^±^)(u,2+- 2 .2 


■» z 

y«^i+y. 


°*(yn+i>-°*(y"r 


n+1 ■'n 


^ “ zAfc 




(5.208) 


.205) 


From (5.208) we observe that; 
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(a) 

(b) 

(c) 


AV„ S 0 

XI 


Since “ ~Z“ both terms in (5.208) cannot vanish 

simultaneously unless = y^^^j = 0 
Since 0(At) 

AVj^ ^ ■' 2 At[(y^)^"H«^(y^)^] + O(At^) as At -+ 0 

The discrete system (5. 199)is Liapunov asymptotically stable. 


Effect of Viscous Damping and Additive Forces 

If to equation (5.168), viscous damping and extex'nal forces are 
added, then using (5. 174) we have 


X + 2zx + 0) X + g(x) - p(t) 


(5.209) 


then if sup jp(t)j = p^, all solutions of (5.209) are xiltimately bounded. 


Let 


V(x,x) - ^ [x^ + 2zxx + U)^x^ + 2z^x^ + 2G(x)] 
V(x,x) > 0 x,x 0 


‘ *2 • 2 * -. 2 * 

V = XX + ZX + 2XX + tu XX + 2z XX + g(x)x 


Using equation (5.209) 

V = -z(x^-HD^x^+xg(x)) + p(t) (x+ zx) 
if xg(x) >0 X 0 

V ^ “ z(x^-Ho^x^) + p^(lx| + z|x|) 

^ - zfx^-to^3i 


Let S be the set x^ + (U 



(5.210) 


(5.211) 


(5.212) 

(5.213) 

(5.214) 
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Outside the set S, V<0 

Let Q be the set V £ e, where e is such that S is a proper subset of 
a. Then, for points (x,:^ outside Q, V< 0 (5.215) 

Starting outside n, V>0, V<0, tlierefore, V decreases and tile 
trajectory must eventually enter Q, and once inside Q, the trajectory 
cannot leave Q since V ^ 0 on an. Starting inside n, V > 0, V is in 
general sign indefinite, therefore V may increase, however, it is clear 
that the trajectory cannot leave Q since Q on 80. Thus all solutions 
of (5. 209) are ultimate bounded in0=(0 + 30). 

Consider now tiie discrete form of equation (5.209) 


-V At ,• 
^n+1 ^n " 2 


Let 


'Si+r'i = - f +“^(yn+i+yn) - (p^+i+p^w 

°^ya+i>-'5(yJ 


- At 


^ii+ryn 


[(yj^+zYn)^ + (z^ + 2G(y^)] 


AV = V , -V 
n n+i 


■ - 2C3(yJ] 

Using equation (5.216) 


n+r 


(5.216) 


(5.217) 


(5.218) 


-S6^ 


AY= -zAt 
n 




2 


+ At 




^n+l-yn 


yn+I+^nV 


Using (5.207) we have; 


AV s; -a At 
n 


(!£^) 


y„j.i +n 




( 5 . 219 ) 


(5.220) 


Let 


h+^ = <f > . !2±^ = <y > 

AV^ s - zAt[<y>^ <y^)^] 


+ 




Let ^S) he the set 


(5.221) 




(5.222) 


Since 

w “ w ~ 0(At) as t~* 0. 

—n+1 — n ^ ' 

Clearly ^S) ^ S as At 0 and the previous arguments can be used to 
show that the solutions of tlie discrete equations (5.216) are ultimately 
bounded, just as the solutions of the continuous equations (5,209). 
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Extension of Energy Conserving Algorithjns to Multidegree -Freedom 
Nonlinear Systems 

Consider the system, of conservative nonlinear differential equa- 
tions 

Mx + Kx + V u(x) = 0 (5 . 223) 

where M, K are N X N symmetric positive definite matrices and u{x) 
is a positive definite potential function. For the system (5.233) 

(5.224) 


We may write equation (5.223) in discrete form as 

y -y = ^ [y +y ] 

— n+1 ~n ^ Uii+i — n-* 


V(x,x) = ^ [x"^Mx + + 2u(x)] > 0*^ 


• • T • T *T / V 

V = X Mx + x Kx +x V^u(x) = 0 


Algorithm A 




(S.2Z5) 


-At 


(yn+r>i^''“(Zn+i)+''!'“(ya)] 


Cross multiplication yields 


= 2- tin in + 2„ + 2 u(5^)] 


(5,226) 


Thus if 
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Vn = i [y^My^ + Zn^Zn ^ ^ 






{5.227} 


Hence the discrete equations (5.225) conserve energy in exactly the same 
way as the continuous time equations (5.2E3). 

Accuracy 

Since energy is conserved, y^,y^, are boiindedfor all n, provided 
Zq'Zq^^® bounded. Thus as 

' Zn+i'Zn'^'^t^*) . 

(^(Zn+l) “’^tZn)) [^^(Zn+l^‘^''^^fZn)^ 

(Zn+l-Zn)'V(Za+l)-'''"<2;n) 




(5.228) 


J 


Hence, as At-* 0, the discrete equation of algorithm A became of 
trapezoidal form, and hence this algorithm is second order accurate as 
At “♦ 0. While this algorithm conserves energy, it has two defects. 

a) It is difficult to use, that is, it is not readily computable. 

b) If At is not small, we have not been able to prove that: 

(Zn+l ® implies that, 

'"(Zn+l) “ “(Zn) 

Thus we are unable to prove that the last term in (5. 225) is bounded when 

y , y , , are bounded. 

•^n^ — n+1 


Algorithm B 


An alternative to the discrete gradient operator in (5.225) is the 
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ope rator 


r 


liU = 


N 

[ 

■l k=l 




IV A 

N L 


^n+l^^n 




where 


h k+1 

yn+l’ ’ 


^"k = "(yi+i'yn+i ■ 

,12 k k+1 k+2 

-“(yn+l-yn+1- • -yn+l-yn ’^a 


(5.229) 


i-1 i i+1 nv 

yn+l'^n ' ' 




(5.230) 


?hus 


in+r^ = ^ ti+i'-ij 
-ij = - ^ 


LiLn+l -in-" 


N 




^^iV(yi+i"^n> 

=^2V(yn+l-yn) 


V 


l^SA^Uj^/(y“^l-y” 




Cross multiplication yields 




N 


i^k 




>- 


(5.231) 


(5.232) 


Now 
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N 


W Z = ’^fen+l) - “ten) 


(5.233) 


i.k 


To see this, consider N = 2 
2 


i 1 -"(y^yn) +"te„+i'yn+i) -“te^.y^+i) 

i,k 

+^(yr,ii>yriT) - 


' n' ^ 


^n+1’ J'n+l^ 




= i [2“te„+i, y4i) - 2”(y^.yn)) 
= “ten+i) -'^te„) 


(5.234) 


Thus, using (5. 234) in (5.232) we have: 


^[yn+l“in+l +2n«^Sn+l + ^nfe^+l)) = (5-23S) 

Thus algorithm. B also conserves energy. If the potential u(y) can be 
expressed in the form u 
N 


u(y) = “i(y^)' 


(5.236) 




where u (r) is a positive monotone increasing ftmction of r, then it 
may be shown that 


N 

1 V . 

W L ^i^ 

k=l 


i ” \ 2 t i\2 

■■n+l) -tej 


SO Vi 


(5.237) 


Application of Algorithm B to System (5.223) with Viscous Damping and 
Additive Forces 


If to equation (5.223), viscous damping and external forces are 
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added, we have 

Mx + Cx + Kx +'V u(x) = p(t) (5.238) 

Then if G is symmetric and positive definite and l^(t)|j is bounded, all 

T 

solutions of (5.238) are liLtimately bounded provided x V u(x) s 0. 

Let 

V = -^ [ss(x'^(C-3M)x) -f (x+zx)'^M(x+2;x)+2u(x)] > 0 

+ x'^Kx ■ (5.239) 

where 2z = smallest eigenvalue of = 0 

V = zx'^(C-zM)x + (x +zx)*^M(x'*+zx) +^V^u(x) +x"^Kx (5.240) 
Using (5 . 238), we have 

V = 2 x"^(C- 2 M)x - (x+2x)"^[(C-zM)x +Kx +V u{x) - p(t)] 

•T *T /X 

+ X Kx + X V u(x) 

t 

- -^’(C-zM)x - z(x^Kx +x"^V ■u(x)) - (x+ zx)*^p{t) (5.241) 

•—* Ji, ■ ““ ““ 

Since 

x^V u(x) s 0 

Xl 

V ^“^(C-zM)x- 2(x"^Kx) + l(x+2x)'^£(t) I (5.242) 

Since 2z is equal to %e smallest eigen value of M ^C, C-zM is sym- 
metric and positive definite, hence the first two terms in (5.242) are 

negative definite. Since p(t) is bounded, and the third term contains 

• • T • T 

X and X linearly, there exists a set S: x Mx +x Kx ^ K such that 

outside of S, V <0. Let Q be the set V ^ C, whei'e C is such that 
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s is a proper subset of a Then outside of fiv <0. and applying the 

arguments used previously^ we see that all solutions of (5.238) are 
ultimately hounded in Q . . 

Now consider the discrete form of equation (5.238) 

^n+r^n = ^fen+l+ij 

^^n+1 ' T ^ 


At 

N 


I ^iV^nii-yn 

k=l 




N 


L 


I ^N^^yn+i-y; 

k=l 


N 

n 


, At , . 

(£n+l+£n) 


J 


(5.243) 


Let 


2 y^Ky^ + 2u(y^)] 

(5.244) 


Then 


^"^n = Vn«-^n 

= ^®ten+r2n)^(C-^M)(jr^^^+y_^) 


+ (2n+l -Vn+^(2nil +2n) '^'^(in+l ty„)) 

*2n+ryn)^K(2^^j+y^) + Z(u(y^^j)-u{y^))] 

Using equations (5.243) and (5.237) 


(5.245) 
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+ At |<yjj+2:y^>'^<Pjj> 1 


where 


<x ) = ■«• [x +x ] 
--n 2 ‘■-’11+1 — ir 


(5.246) 


(5.247) 


Since ~—n At "♦ 0, it is clear that the right hand side of 

(5. 246) tends to 

-At [y'^(C- 2 M)y + zy*^Ky ] + At l(y +zy )*^p 1 as At -♦ 0 

(5,248) 

Thus, for small At, the con-fcinuous time system (5,238) and the discrete 
time system (5. 243) hehave in essentially the same way and the solutions 
of equations (5.243) are ultimately bounded. 

Accuracy 

Using the fact that the solutions of (5. 244) are ultimately bounded, 

it is easily shown that the discrete gradient operator used in (5. 244) has 

the following form as At 0 
N ' 

1 1 


1 

k=l 




N 

k=l 


I =i[Vu(y^^^)fVa(y^)3 




+ OiAt^^) 


(5,249) 




Thus, as At"*0, equations (5. 244) are of trapezoidal form, there- 
fore algorithm B is second order accurate as At tends to zero. 

Algorithm B conserves energy and has some nice stability 
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properties ; it -anfortunately has two major defects 

a) It is difficult to use, that is, it is not readily computable. 

b) The property (5.237), which is necessary in order to prove 
the ultimate boundedness properties of equations (5.243), is valid only 
for a restricted class of potential functions u(y). 

For this reason, we now turn to an alternate formulation using 
Lagrange multipliers . This formulation was suggested by my colleague 
Dr. T. J. R. Hughes and was developed jointly with him and his student 
Mr. W.K. Liu. (3) 


Algorithm C 

Consider the system of conservative nonlinear differential 
equations 

Mx + Kx + V u(x) = 0 . (5.250) 

where M, K are N x N symmetric positive definite matrices and u(x) 
is a positive definite potential function. We know that (5.250) has the 
Liapunov function 


V(x, x) = ^ [x^Mx -f x”^Kx + 2u(x)] = constant 

If we write (5.250) in trapezoidal discrete form 

At r * • T b 

^(in+1 ^ +5^) + +7U J ^ 

We now wish to constrain (5.252) such that 

1 • T • IT 

-Pj y ,.My ,i-i--^y ,Ky ,+u , 

2 ^n+1 An+1 2 Ln+1 ^ n+1 


(5.251) 


(5.252) 


= const 


(5.253) 
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Using the first of equations {5.252) 

2 • 

^ fln+1 "I.J ” In (5.254) 

Substituting into tlie second of equations (5.252) 

2 

= “ (^) (5.255) 

which may be written 

+ (f - (f )"[Ky^+7u^] 

(5. 256) 

Using (5.254), equation (5.253) maybe written 

= Z fexx.i -Zn- f -y^- ^ i) 

+ (^) [z2n«K2n+l+’'n+i] 

■ [zI^Ky^+uJ = 0 (5.257) 

het us now construct the functional 




■n+r'Zn+i+^n+1 


) 


■^2„+iB<yn+V’iJ {^f (S.ZSS) 


Then 


6y, 


T 8F 




— n+1 8y 

— n+1 


Myn+1+ V [K2n«+’“n+ll 
2 


- M[y„+At_J^] + (^) tKy^+Vuj] 


(5-259) 
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Thus necessary and sufficient conditions that equation (5.255) hold, are 
that (5.259) vanish for arbitrary ’’variations” • 

In order to force equations (5.252) to conserve energy we com- 
bine (5.250) and (5.257) tlirough the use of a Lagrange multiplier our 
new functional is 




(5.260) 


Thus 


T r 

[ 




^2n+l 


®2n+l 




(5.261) 


If (5.26l) is to vanish for arbitrary 6X then equation (5.257) is 

satisfied and in addition; 

"^2n+l + (^) ) [Ky^^i+VuJ 

.2 




Thus 


^.2 


(l+X)[My^^l+ (^) [Ky^^j+Vu^^j)j 

= (l+>.)[MyJ + At (l + I) My^ - (^)^(Ky^tVu^) (5.263) 

The new algorithm thus consists of the two systems of equations, (5, 26l) 
and (5. 26,3). 

» 

When y has been determined, y can be calculated from 
— n+1 — n+1 

equation (5. 254). 

The new algoritlrm is solved using a variant of the Newton- 
Raphson method. 


If i denotes the iteration number tlien if we define 
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„i+l „ ^ _ A i 

-Ji+l “ ■ ^2n+i > 


= AX* 


We have 


(X^l) [MAy^j + 

+AxfM(2^-yJ + (v)Vy4, ^ Mi J 

*■ (^) - (i+^)y„At = 0 


i T 

^Za+> 






(5.264) 


(5.265) 


+' 3 (z 4 l) = 0 


where 


(5.266) 



9y \ ' 

— n+1 


Equations (5.266), (5.267) are of the form 


= b" 


where 



(5.267) 


(5.268) 


(5.269) 


where 
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A^2 = ["<4l- ^ - f^n) " ''"n+l” 


4l = <42>" 


+ (^)^ T% 


^ “ <^(^+7 ) 

Rewriting equation (5.268), 

^21 ^4+1 "" *^2 

+ a] 2 AA^ = b^ 

From the second equation in (5.271), 

“ ^^11^ i ” ^12 

Using the first equation in (5.271), 


4l ^^+1 y " 4l^^ll^ ^ ^12 *^2 

AA' = [{Ai2)Vi7)“'' bi - bj]/(Ai25^<^]\)“’' A^p 


= U(\^yhh\ - aa'a]^] 






(5.270) 


(5.271) 


(5,272) 


(5.273) 


(5.274) 


i i 

It should be observed that and AA can be obtained with only one fac- 
torization of a|i and two forward reductions/back substitutiors. Thus one 
additional forward reduction /back substitution is required when compared 
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» 


% 


with the Newton- Raphson implementation of the trapezoidal algorithm. Thus, 
unlike algorithms A and B, algorithm C is readily computable. 


Accuracy 

Using equations (5,252) and (5.263), the new algorithm may be rewrit- 
ten as: 


^+1 " ^ 2 ^-^n+l 

inn-T in = - ? ^ ’“n+1 * 

f -y 

? ^+1 '^^+1 * r Ij ^ n+l '^+1 




) (5.275) 






From the first two equations 

J [y^+i Min^ ^ r 

+ ? <Vl - <%+l - Jt 

(5.276) 

Using the third equation of (5.275) 


m '2n+l - in’’’ ■ It 

= T (Vl - (VUn^.1 + VUn) + u„ - Un+, - (5.277) 


Thus 
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A 

1+A 




where 

'N 

c = ^[J(Si) - a(E2)] 

0 < 0^ < 1 , i = 1 ,2 

Hence, equations (5.275) may be rewritten 

^+1 ~ ^ T ^^+1 

in^-l - in = " ¥ ^ ^ / 


(5.278) 


(5.279) 


(5.280) 


If u has continuous first and second partial derivatives, then|!C||^ 
is bounded and we may apply standard techniques to (5.280) and show that; 

116^11 = llin • ^11 11^ ■ 2r.ll- “ ‘ as At-0 (5.281) 

Thus algorithm C is also second order accurate as At 0. 

The Lagrange multiplier technique is clearly superior to the other 
techniques, and while in the present analysis the constraint was that of 
conservation of energy, the technique can be used with any appropriate 
constraint. For example, if the technique is applied to equation (5.238) 
and we wish to ensure that the solutions of the discrete equations will 
be ultimately bounded, given that the solutions of (5.238) are ultimately 




V 


I 
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bounded, the appropriate constraint would be the discrete form of (5.241) 
obtained by integrating from t^^ to t^^ + At, i.e., 

Vl ■ ''n = - X + iJ(G - zM)i„ 

+ z(4l 

■<■ (Vl ■*■ zVl >^ fin+1 'in 

Replacing by (y^^_j_^ - y^) ■* ^ , the constraint equation becomes 

^ '^Vl - - in>^ M(4(VT^’ - in’ 


- ^ [ziJ(C -zM)i„ 'in '^'in ^ 

+ K^n + 2u(j!^)] 

f "n 'Vl - in’^ (G-2H)(4 'Vl - in’ 
■ + 4(C - zM)y„ + l<Vi + 1%) 


* 2(4+1 ''%+! + 4 %’ 



' ®'Vl’ = ° 


(5.283) 


The functional oi" equation (5.253) is replaced by: 


= ^4+1 %-H 'f) yJ+1 %-f-i 
+ (t>^ L 5- 4^-1 >^+1 ^ Vl^ 


- 4+1 ■<• it in) + f 

+ (f)^ 4+1 t 7u„^, - p„^., - p„] 


(5.284) 


The variation of the functional (F(y^^^) + ^G(y^^-[)) yields the new algo- 
rithm: 


9F ^ , 3G 

»i^+1 ®J!n+1 


0 


G(i„+l) = 0 


^ n +1 = h (^+1 ■ ■ 4 


(5.285) 


It may easily be shown that (5.285) is second order accurate as 'At ■»■ 0. 
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6, Application to the D.ynamicaT Analysis of Large Space Vehicles 

Consider the system of differential equations obtained by applying 
finite element techniques (or other techniques) to some complex space 
vehicle- The equations are likely to be of the form: 

M^+ Cx + K>c + v^u = £{t) (6.1) 

M is an symmetric positive definite matrixj C and K are N^N sym- 
metric positive semi-definite matrices, and u(^) is a positive semi-definite 
potential function. (In the present analysis we shall neglect terms in 
(6.1) which arise due to steady rotation of the vehicle*) 

Since one of the primary objectives of any structural analysis is to 
determine the stresses in the vehicle, it is desirable to make N, the number 
of coordinates, as large as possible so that stresses may be determined 
accurately. The number, R, of modes of the structure exhibiting signifi- 
cant response is usually much smaller than N. This poses a serious diffi- 
culty for the direct numerical integration of (6.1), since, as we know, 
the accuracy of any numerical scheme is determined not by the time step 
At, but by {(o.At) where w. is the highest “frequency" which can be excited. 

J U 

In order that the higher modes be integrated accurately. At may have to be 
very small, much smaller than is either practical or economically feasible. 
In the case of linear systems, the use of algorithmic damping or post- 
filtering successfully overcomes this difficulty by suppressing the higher 
modes, which are inaccurately integrated when a moderately small value of 
At is used, thus resulting in reasonably' accurate representation of the 
lower modes. 
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Hany analysts use algorithmic damping or post-fi 1 taring to achieve 
the same result for nonlinear systems. Care must be exercised when doing 
this, since nonlinear systems can exhibit internal resonance, a phenomenon 
in which higher modes, though not excited by external forces, can be excited 
by nonlinear coupling to the lower modes. To illustrate this phenomenon, 
consider the following problem: 


x^ + 2z^x^ + + uCx-j-Xg)^ = P-jCOS cjjt + PgSin ut + P^cosSmt 


> ( 6 - 2 ) 


Xg + 222X2 +9x2 + u(x2~x-j)'^ - 0 


J 


where 

p2 = “2z^wA 
P3 = I vP? 


( 6 . 3 ) 


Since the second equation is not excited externally, it may seem reasonable 
to set X 2 = 0; the fi'rst equation then has the solution 


C" 


4 


x-| = A cos wt (6,4) 

If we now turn to the second equation, regarding Xg as small, 

Xg + 222^2 ^ ^2 " ^^1 “ "T 4 ^^cos Suit (6.5) 

If (D -V 0(1), the first term on the RHS of (6.5) causes no trouble, 
however, the second term will cause resonance, ad if Z 2 is small, will cause 
significant response in X 2 » We therefore see that even though the second 
equation in (6.2) is not externally excited, it car. still be driven by the 
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first coordinate This is a simple example of internal resonance. 

Quite frequently as a preliminary to performing a nonlinear dynamic 
analysis, a modal analysis of the linearized system will be carried out. 

The modal analysis can be a very useful tool in structuring the nonlinear 
problem for dynamical analysis. 

As a preliminary, let us first put equations (6.1) into canonical form 


let 

Let 


y X 

C = 

JC 

U(x) ^ V(y;) 

= £(t) 


Using (6,6) and (6.7) in (6.1) 




ly + +JC^ + VyV(y) - £(t) 


y(o) = A i(o) = b 

Let T be the orthogonal matrix which di agonal i zes JC. 


( 6 . 6 ) 


(6.7) 


( 6 . 8 ) 


Let i “ T ^ 

tItct = A 

V(i) = W(z) (6.9) 

m = 

f(t) = (f,(t)} 
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Unless (6.1} has classical normal modes, S)is not diagonal. Using 
(6.9) in (6.8), 


I ^ + S)z + A ^ W(^) = f^(t) 


z(0) = a 


z = b 


Suppose that lf-(t)[ < e for j > P , e « 1. 

J 




If we suppose that the j ^ modes, j > P are at most weakly excited, let 
us set 2 . = 0, j C(P+1»N). Then, 

J 


W(z) = W(^) 
equation (6.10) becomes 


I il = fq(t) 

I ^ -^2-2 ^2^^^ 


. V represents the rigid body modes 


-2 


represents the first (P-6) flexible modes. 


T 

^ 22 22 X (P-6 ) clamping iTietric associated with the 


^ modes. 


V 


r 2 




(07 


7 


/ 


h 

■> • ; 

♦ 


“i 

% 

4 

> 

K 

■ • 

*2 

0)n 


P 


P 







(6.T4) 


If P is not too large, we can select a At such that (uipAt) v o.l , 
then using any of the "energy conserving" altorithms of Section 5, equations 
(6.13) can be integrated with good accuracy. Having determined one can 
easily compute the pliysicaT coordinates, 




= 


h 


where T is the NxP matrix having as its columns the first P eigenvectors 
of the linearized problem. To check if there is any significant response 
in the neglected modes, due perhaps to internal resonance, we approximate 
the remaining modes by the system of uncoupled equations 








+ 



= 

3z. 



(6.15) 


j e(P-M ,N) 


We note in passing that equations (6.15) will be exact if the linearized 
part of (6.1) has classical normal modes. Not all the modes in (tl.15) need 
be examined’, only those for which 



« 1 


(6.16) 


where p and q are integers and Aj^ is an element of Ag. If no mode of (6.15) 
shows significant behavior, we can be reasonably sure that the solution of 
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equations (6.13) will give a reasonably accurate representation of the 
solution to equations (6.1). 

If any mode of (6.15) shows significant behavior, we can be reason- 
ably sure that the solution of equations (6.13) will not give an accurate 
representation of the solution to (6.1). In this case, the modes of 
(6.15) which show significant behavior must be included in the solution 
of the problem. This presents a serious problem in the general case, 
since we require that (wj^At) ^ 0.1 for accurate integration of the system. 
If only a few modes of (6.15) show significant behavior, it may be pos- 
sible to treat the problem in an efficient manner. 

"k 

Let = 

where k, m are the modes of (6.15) showing significant behavior. 



Let 



(6,18) 


Equation (6.10) may, in this case, be written 




I +^^^22^2 '*'^23^3 ^2^ ~ 

I h * ® 23^2 ®33j-3 ''32.3 = ° 


(6.19) 


J 


This first set of equations is integrated using a At^ appropriate to the 
highest eigenvalue in Ag. The second set of equations is integrated using 
a Atg appropriate to the highest eigenvalue in A3, say Atg = At.| , K an 


integer. The values of ^ appearing in the second set of equations can 
be obtained by interpolation from the solutions of the first set of equa- 
tions. 

Internal resonance occurs most often in systems where the eigenvalues 
of the linearized system are Integrally related, and where the nonlinear 
system is subjected to a steady state single frequency excitation; for- 
tunately these two situations do not appear to arise too frequently in 
the space vehicle problem. Nevertheless, such situations can arise, and 
the analyst should be aware of them. 



Appendix 1 


Generalization of Theorem 8 


Theorem . Given the linear difference equation 

lA(n)| / 0, ||A(n)!l <« 

then A1 is uniformly Liapunov asymptotically stable at x = 0 iff there 
exists a bounded, symmetric, positive definite matrix P(n) such that, 

T 

1) P(n) = P (n) positive definite and bounded above & below 

ii) A”*^(ii) P{n+1) A(n) - P(n) f - 0(n) 

T (A2) 

iii) e(n) = 0 (n) positive definite 

iv) lle(n)ll < M 2 (n^) V' n > n^^ and ^ 


Proof Sufficiency 

Suppose that there exists such a matrix P(n) satisfying A2 

Let V = P(n)x (A3) 

n — n — n 

Since P(n) is positive definite and bounded. 


Vi = 4+1 

Using (Al), 

= xj A'^(n) P(n+1) A(n)x,^ 

' <Vl ■ ''n> ' 


(A4) 


(A5) 

(A6) 
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Using (A2) 


AV„ = - 4 e(n)x„ < 0 


(A7) 


1 


I 


f 




'n+1 


< < V 


n-1 




(A8) 


Since P(n) is bounded Vn , is finite if H^^l! is, and since is 
zero only ii^ = 0» lienee and therefore ii^li tends to zero as n 

tends to infinity. Since the result is independent of n^, the trivial solu- 
tion X = 0 is therefore uniformly Liapunov asymptotically stable. 


Necessity As in the proof of necessity for Theorem 8, it is easily shown 
that P(n) satisfies equation (4.53), thus: 

P(n) = I *(j,n)^ 0(a) $(j,n) (A9) 

j=n 

Thus if (A1 ) is uniformly asymptotically stable at the origin 


il$(j,n){UM^ ^(j-n) 

0<6<1 , j>n, \/n>n^ 

Using (A2iv) 

1) ||P{n)|| < j m 2 «2U-n) 


1 - 5 " 

11) Phn) - (I fhj.n) e(j) ♦(j.n))'^ = P{n) 
j'n 


(AlO) 


(All) 


(A12) 

(A13) 
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iii) x^P(n)x = I (^(j>n)x)^ e(j}(^(j,n)x) >0 

j=n 

if '^(j»n)_x r 0 

and since |#(jjn)l 5 ^ 0 , 

P(n)^ > 0 X. 7 ^ 0 (A14) 

Thus If (Al) is uniformly L.A.S. at 0, there exists a matrix P(n), 
symmetric, positive definite, and bounded which satisfies (A2ii). 

iNote : It is clear that P(n) is difficult to compute, except through the 

use of {A 9 ), which requires the unknown transition matrix ®(j,n). Its main 
use is in proving theorems. 
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